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ABSTRACT 

Acoustic propagation in an atmosphere with a specific form of a temperature 
profile has been investigated by analytical means. The temperature profile used is 
representative of an actual atmospheric profile and contains three free parameters. 

Both lapse and inversion cases have been considered. Although ray solution have 
been considered the primary emphasis has been on solutions of the acoustic wave 
equation with point source where the sound speed varies with height above the ground 
corresponding to the assumed temperature profile. The method used to obtain the 
solution of the wave equation is based on Hankel transformation of the wave equation, 
approximate solution of the transformed equation for wavelength small compared to the 
scale of the temperature (or sound speed) profile, and approximate or numerical 
inversion of the Hankel transformed solution. The solution display the characteristics 
found in experimental data but extensive comparison between the models and 
experimental data has not been carried out. 
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1.0 ltjl PO DUC TIQ N 

Although the propagation of acoustic signals through the atmosphere has been 
studied for many years, and most atmospheric effects are understood in the qualitative 
sense, quantitative modeling of most of these effects has become an area of interest 
only recently. The dominant effect occurring in atmospheric propagation is the 
spreading of acoustic energy associated with a wave propagating in three dimensions 
over an ever increasing area, the well known spherical spreading effect, which occurs 
in an isothermal, unbounded atmosphere. In addition to this, acoustic waves are 
absorbed by the atmosphere, reflected and absorbed at the ground surface, scattered 
by turbulence, and refracted by both wind and temperature gradients. 

This report summarizes a project to develop models for the propagation of 
acoustic signals from a point source above a finite impedance ground surface in the 
presence of temperature gradients in the atmosphere. The situation of interest is the 
case of sound from a source located within a few meters of the ground propagating to a 
receiver located within a few meters of the ground through the temperature gradient 
that commonly occurs just above the ground surface. Best[1], Gieger[2], and 
Reynolds[3] all discuss the temperature gradient in this region. Within one to two 
meters of the ground the temperature generally goes through a diurnal cycle with a 
lapse condition, temperature decreasing with height, occurring in the afternoon and an 
inversion condition, temperature increasing with height, at night, see Figure 1.1. 

Shortly after sunrise and sunset the atmosphere goes through a nearly isothermal 
period when the transition from lapse to inversion or inversion to lapse condition is 
under way. This simple picture of the very complex atmospheric dynamics near the 
ground can be upset by significant winds which increase the mixing near the ground 
surface and tend to lead to a more isothermal situation, or to an overcast which can 
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prevent a strong lapse condition from developing by blocking the insolation or prevent 
an inversion from occurring by blocking the radiation from the ground to the night sky. 

References 1 ,2 and 3 all discuss the classic logarithmic temperature profile given 
by 


T 



(1.1) 

which is based on empirical results. This result, although fitting the experimental data 
well, certainly is not reasonable either for heights very near the ground or very far 
above the ground. In addition, logarithmic functions are generally more difficult to deal 
with in an analysis then are algebraic functions. For this latter reason the profile used 
in this study is 


T = T — 5L_ 

00 1 + az 

(1.2) 

This form of the temperature profile is shown in Figure 1.2 along with some 
temperature data obtained by Butte rworth[4j. The agreement between the data and 
the assumed function fitted to this data is excellent. Also as compared to (1 .1 ) the 

physical meaning of the parameters in (1 .2), T«, AT, and a, are clear. The assumed 

temperature profile asymptotically approaches the temperature T^ high above the 

ground. At the ground the temperature is T^ + AT, thus the change in temperature 
between the ground and far about the ground is AT. The derivative of temperature 
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with respect to height evaluated at the ground surface is - a AT. Thus 1/a is the scale 

over which the temperature change AT occurs. For example at a height z = 1/a 

one-half of the total temperature change AT has occurred. The temperature profile of 

(1 .2) can be used to represent either a lapse or inversion condition. For a lapse 

condition the parameter AT is positive and for an inversion it is negative. 

The equation governing the wave motion is the simple acoustic wave equation 
with a sound speed varying with height and with a point source term, 

_L_i£.V 2 p = i e i "'5(r)5( 2 -s) 
a(z) dt n 

(1.3) 

At the ground surface, z = 0, a normal impedance boundary condition 


P-T^t 

ipco ° z 


(1.4) 


is assumed. High above the ground, z -> «, only outgoing waves are permitted, a 
radiation condition. At the source height, z = s, the pressure field is to be continuous, 
and to satisfy the conditions implied by (1 .3). 

Section 2 contains a discussion of the acoustic rays that characterize both the 
lapse and inversion cases for the assumed temperature gradient, (1 .2). This both 
yields a quantitative understanding of the propagation phenomena, and plays an 
integral part in understanding the modeling that follows. The model for the lapse 
condition is developed in Section 3 and the inversion model is described in Section 4. 
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Section 5 contains a discussion of the conclusions developed during the project. 
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2.0 ACOUSTIC RAYS 


Acoustic ray tracing is a relatively simple procedure for an axisymmetric case 
which yields a great deal of qualitative information about a given propagation situation. 
Only a brief discussion is given here. More detail is given in [5]. For an horizontally 
stratified atmosphere the acoustic rays may be determined from an acoustic form of 
Snell's law 


Cos 8(z) Cos 6(s) 
a(z) " a(s) 


( 2 . 1 ) 

where the source is located at the height s. Here 0(z) is the angle between a ray and 
the horizontal at the height z, see Figure 2.1. Thus the right hand side of (2.1) is a 

constant for a ray emitted from the source at an initial angle 0(s). Using (1 .2) to obtain 


a 2 (z) = a 2 ( 1 + 

oo 


AT 1 
T „(1+az) 


) 


( 2 . 2 ) 

which describes the sound speed as a function of height, the slope of a ray initially 


emitted from the source at an angle 0(s) is determined from (2.1) to be given by 


dz / Az + B* 
dr “ ± \J Cz + D 

(2.3) 
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where 


A = a[ 1 + as + Y^-(1 + as) Cos 2 0(s) ] 


B = [ 1 + as + Y“* (1 + as ) ( 1 + y" ) Cos 2 0(s) ] 


(2.4) 


(2.5) 


C = a ( 1 + as ) Cos 0(s) 


(2.6) 


and 


D = ( 1 + y~) ( 1 + as ) Cos 2 0(s) 

oe 

(2.7) 

Equation (2.3) can be integrated to obtain the ray paths. Different results are obtained 
in the lapse and inversion cases and these will be considered separately in the next 
two sections. 

2.1 Lapse Case 

In the lapse case the quantity A is positive and integrating (2.3) to obtain the rays 
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yields four different cases. For rays going upward from the source (using the positive 
sign in (2.3)) 


r = F(z, e(s))-F(s. 0(s)) 

for rays going downward initially from the source (using the negative sign) 

r = F(s, 0(s))-F(z, 0(s)) 


(2.6) 


(2.7) 


for rays that were initially going downward but have been reflected upward at the 
ground (using the positive sign and (2.7)) 

r = F(z, 0(s)) + F(s, 0(s)) - 2 F(0, 0(s)) 

(2.8) 

and for rays that initially were going downward and were refracted upward before 
reaching the ground (again using the positive sign and (2.7)) 

r = F(z, 0(s)) + F(s, 0(s)) 


(2.9) 


The function F(z, 0(s)) is given by 




( 2 . 10 ) 
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where 


and 


E = A D - B C 

= a ( 1 + a s + ) ( ~ ) ( 1 + a s ) Cos 2 0(s) 

— M 

( 2 . 11 ) 



C ( A z + B ) ' 
| A | (Cz + D) 


( 2 . 12 ) 

The absolute value of A in (2.12) is immaterial in the lapse case where A is always 
positive but significant in the case of an inversion where A may change sign. 

The rays are identified by the parameter 6(s), the initial angle at which the ray 
leaves the source. Thus a ray initially propagating downward and identified by a 


particular value of 6(s) will either be reflected upward at the ground or refracted upward 
at a turning point. In either case the reflected or reflected ray will be identified by the 


same value of 0(s) as the initial ray. Figure 2.2 is an example of the rays calculated 
from (2.6) to (2.9). 

Setting 0(z) = 0 in (2.3) yields an expression for the height at which an initially 
downward propagating ray becomes horizontal, the turning point, as 
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(2.13) 


B 

Z »P a! ' A 

Solving this expression for Cos 6(s) yields 


Cos 6(s) = 

(2.14) 

which identifies the ray having a turning point at a height z^. The ray that grazes the 
ground and is the boundary between reflected and refracted rays can be found by 
setting = 0 in (2.14). The ray that divides the initially upward and downward 

propagating rays is identified by 0(s) * 0. 

For the ray that grazes the ground, and is identified by the value of 6(s) defined by 
(2.14) with ZipsO, the function F(0, 0(s)) = 0 and thus this ray is defined by either (2.8) 
or (2.9). Similarly for 6(s) - 0 the ray can either be obtained from (2.6) or (2.9) since on 

this ray F(s, 0) = 0. At a turning point F(Ztp, 0(s)) = 0 when 0(s) is given by (2.14). 

The shadow boundary is more difficult to locate. It is bounded by refracted rays 
(2.9) and at a fixed height z is the maximum possible value of r for rays with turning 
points below that height. Thus the ray tangent to the shadow boundary or caustic at the 
height z is identified by solving the equation 
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0 


dr _ aF(z,e(s)) aF(s, e(s)) 

9 0(s) ~ ae(s) + ae(s) 

(2.15) 

for 0(s) and then using that value in (2.9) to determine the location. The derivative in 

(2.1 5) can be obtained from 


if. 

dy 


+ 


{(1 + az)(1 


o (1 - y 2 ) 2 
+2Y 2 ] 

OO 

i 

2[ 1 -T 2 ] 2 


JZ'a 




2 , „ AT . 

z-y (1 + az+y-) 






where 


(2.16) 



Cos 6(s) 


(2.17) 

Due to the complexity of this expression an analytic solution is not possible and either 
a numerical solution of (2.15) must be obtained or the approximate relations given in 
[5] must be used. 
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2.2 Inversion Case 


As in the lapse case several different types of rays occur in the inversion case. 

Integrating (2.3) with AT '/T m negative also yields three different forms for the function 

which determines the rays depending upon whether the quantity A given in (2.4) is 
positive, negative or zero. The meaning of these three cases is discussed below. For 
rays that are initially angled upward (using the positive sign in (2.3)) 


r = F.(z, 0(8)) - F.(s, 0(s)) 


(2.18) 

where i = 1 , 2, or 3 depending on the initial angle of the ray leaving the source. Rays 
with i =1 leave the source a sufficiently large angle upward so they do not have a 
turning point and are never refracted downward. The case i =2 corresponds to the 
limiting ray that has a turning point at infinite height. Rays described with i = 3 have 
turning points at finite height and are alternately refracted downward and reflected 
upward at the ground. These are the rays trapped by the inversion. 

Rays that initially are angled downward are given by (using the negative sign in 

(2.3)) 


r=F.(s, 0(s))-F j (z, 0(8)) 


(2.19) 

for all three cases before they are reflected upward at the ground. The reflected 
waves are given by (using the positive sign) 
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r = Fj(z, 0(3)) + F.(s, 0(s)) - 2 F.p, 0(8)) 


( 2 . 20 ) 

in all three cases. Note that after reflection the i = 3 rays are refracted downward and 
reflected upward from the ground repeatedly. In the case of these i = 3 type rays four 
more forms exist. For rays that initially were angled upward (using the negative sign) 


r = -F 3 (z, 0(s)) - F 3 (s, 0(s)) - 2 n F 3 (0, 0(s)). 

( 2 . 21 ) 

after they have been refracted downward and have been reflected n times from the 
ground. For rays that initially were angled upward and have been reflected upward n 
times from the ground and have not been refracted downward following that reflection 


r = F 3 (z, 0(8)) - F 3 (s, 0(8)) *2 n F 3 (0, 0(8)) 

(2.22) 

Thus an i = 3 type ray leaving the source upward is first described by (2.18) or (2.22) 
with n = 0 before it is refracted downward through a turning point. After it is refracted 
downward the first time it is described by (2.21) with n = 0. Following its first reflection 
from the ground it is given by (2.22) with n = 1 , then by (2.21) with n =1 between 
refraction through a turning point and reflection, then (2.22) with n =2, etc. 

For rays that initially were angled downward (using the positive sign) 


r = F 3 (z, 0(s)) + F 3 (s, 0(3)) -2 n F 3 (0, 0(8)) 


(2.23) 

after they have been reflected upward n times from the ground. For rays that were 
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initially angled downward (using the negative sign) 


r = - F 3(2, 0(s)) ♦ F 3 (s, 0(s)) -2 n F 3 (0, 0(s)) 

(2.24) 

after they have been reflected n times from the ground and have been refracted 
downward through a turning point. 

Thus an i = 3 ray that is initially angled downward at the source will first be 
described by (2.19) or (2.24) with n = 0 until it reflects from the ground, then by (2.20) or 
(2.23) with n = 1 between reflection and refraction through a turning point. Following 
the turning point and before the second reflection (2.24) with n = 1. Then by (2.23) with 
n = 2, etc. 

The functions Fj are given by 


^(z, 0(8)) 


J (Az + B)(Cz + D) E 


In 


2A/AC X 4 




(2.25) 


for i = 1 , 


F,(z, 0(8)) = 


3C/B 


(Cz+D) 


for i = 2 , and 


(2.26) 
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F 3 (z,e (s ,) J( Az + B >( £*+£> , 


a/ 


==Tan' (-) 
AC 


(2.27) 

for i = 3. As discussed above, the i = 1 case occurs for A greater then zero or for 
Cos 0(s) < Cos 0 (0 > 0) where 


| 

| 



( 2 . 28 ) 


Rays with values of the initial angle, 0(s), greater the value of the angle given by (2.28) 
then escape from the trapping effect of the inversion. The ray with A equal to zero or 

0(s) = © is the limiting ray that has its turning point at infinity (the i = 2 case), while the 

rays with i =3 correspond to negative values of A or Cos 0(s) > Cos0 (0(s) < 0), and 
are the rays trapped by the ground. Rays with initial downward slopes can be divided 
in a similar manner but in all cases at least one reflection occurs before the ray 
escapes the trapping effect of the inversion, is the limiting case, or becomes trapped by 
the inversion. 

Figure 2.3 is an example of the rays calculated from the above equations for the 
case of an inversion. 
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3.0 LAPSE CASE SOLUTION 


The solution of the problem posed by equations (1 .2) to (1 .4) in case of AT > 0, a 
lapse condition, was undertaken first. The general approach used in both the lapse 
and inversion cases was to first separated out a sinusoidal time dependence from the 
pressure, and then to Hankel transform the governing equations with respect to the 
horizontal distance from the source, r, to reduce the number of independent variables 
to one, the vertical height, z. This reduces the governing equation for the transformed 
independent variable to an ordinary differential equation for which an approximate 

solution can be obtained. This solution contains the Hankel transform variable, |3, 
which replaced the horizontal distance in the transformed governing equation. The 
transformed solution must then be inverse transformed to return to physical space. 
Because of the complexity of the solution this inverse transform can not be carried out 
exactly and either an approximate inversion must be used or the inversion must be 
carried out numerically. Both approaches were used in the lapse case, in both of 
these approaches it is necessary to to interpret the Hankel transform variable as a 
complex variable and to continue the solutions off the real axis for the transform 
variable. This is not an intuitive process as the physical interpretation of the transform 
variable is lost off the real axis. This process will be discussed in detail below. 

With the solution obtained for complex values of the transform variable attention 
will be turned to the inversion of the transformed solution. The methods used are the 
classical saddle point approach and a Fast Fourier Transform (FFT) based numerical 
method. These methods are describe in detail elsewhere and will be only described 
briefly here. Finally the results of these approaches will be described. 
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3.1 Transformation and approximate solution 


The time dependence in the governing equation and boundary conditions, 
equation (1 .3) and (1 .4) can be removed by assuming 


p(z, r, t) = e** 5(z, r) 


(3.1) 

The Hankel transform or two-dimension Fourier transform for an axisymmetric function 
can be defined [6] as 


G(z,P) = J G(z, f) r J 0 (Pr) dr 
0 


and the inverse transform as 


(3.2) 


G(z, r) = J G(z, p) |3 J 0 (fJr) dp 
0 

(3.3) 

The use of transform methods in solving partial differential equations arises from the 
fact that an appropriate transform will convert a particular type of derivatives into an 

algebraic term expressed in terms of the transform variable (p in (3.2) and (3.3)) in 
place of the original physical independent variable (r in (3.2) and (3.3)). Thus the 
number of independent variables in the partial differential equation will be reduced by 
one and the transform variable acts only as a parameter in the transformed solution. In 
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the case of Hankel transforms the radial'aependence in the Laplacian operator 
expressed in cylindrical coordinates for a axisymmetric function is converted to an 
algebraic term, see [6] for more details. Applying (3.2) to (1 .3) leads to 


. 2 l 2 
dz a 


P 2 ] G = - — 5(z-s) 

AT 2k 


(3.4) 

The term on the right hand side represents the source. The homogeneous form of this 
equation would have a solution with an oscillating behavior if the term in square 
brackets was positive and an exponential behavior if it was negative. Thus the 
transition occurs when the term is zero or 



Notice that if 


(3.5) 



Cos 6(s) 


(3.6) 

then equations (3.5) and (2.14) are very similar and we can interpret (3.5) as giving the 
value of P that causes the turning point (the location of the transition from oscillatory to 
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exponential behavior as well as a horizontal ray) to be located at the height z. But (3.6) 
remains to be be interpreted. From (2.1) and (2.2), however, one finds that p is equal to 
(c a/a„ ) times the cosine of the angle that a ray, which was initially at an angle 0(s), 
makes in the limit as height tends to infinity. Thus just as we have used 0(s) to identify 

a ray we can also use the Hankel transform variable p. This discussion emphasizes 
the close relationship between the model being developed and the ray description. 
These parallels will be also be pointed out below. 

If height is nondimensionalized in (3.4) using the scale of the temperature 

gradient, a, a uniformly valid approximate solution to the resulting equation can be 

obtained, using the method presented by Nayfah [7], for large values of <a/(a„a) as 


G = ■ y - A - . h 1 (n(z, P)) + E L..- h 2 (T|(z, p)) 

>/ 9 2 (z. P) yj g z {z, p) 

(3.7) 

where 


3 

g 2 (z, p) = (i +oz + 


AT 

T 
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(3.8) 


and 



3 

n(z.P) = (|0 g(z.P) 


(3.9) 


(3.10) 


g 2 (z. P) = 


3g(z, P) 

dz 


(3.11) 


The modified Hankel functions h., and h 2 are defined in [8] by 
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h,©=(§) ¥ hV'( 1 5 2 ) 

3 

(3.12) 
and 

h 2 ©-(|) 5 ? h®(|5») 

3 

(3.13) 

A and B are constants that have to be determined. This solution is rather complex and 
is expressed in terms of unfamiliar functions. Several important features of this solution 
need to be discussed to understand it. 

The functions h,(^) and h 2 (£) have a complicated behavior [8]. For real values of 
the argument both h 1 and h 2 yield a complex results which is oscillatory with an 
algebraic decay of the amplitude for increasing magnitude of the argument. The 
function can be shown to represent downward propagating waves (for e it0 * as used 
here) and h 2 upward propagating waves. The oscillatory behavior also occurs for hj 

when the phase of the argument is equal to 2 tc/ 3 and for h 2 when the phase of the 

argument is -2n/3. When the phase of the argument is n/3, h, decays exponentially 

and h 2 grows exponentially. When the phase is -n/3, h 2 decays exponentially and h 1 
grows exponentially. 

The solution of the point source problem is closely related to the plane wave 
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problem discussed in [9]. In wave problems where transforms have been used the 
inverse transform can generally be interpreted as a superposition of plane waves over 
a range of angles. In the case of the Hankel transforms used here the limits of 

integration on p are from zero to infinity and p can be interpreted as co/a,,. Cos 0„ where 

0„ is the angle between a ray and the horizontal in the limit of height tending to infinity 

as given in (3.6). Thus a value of p of zero corresponds to a ray which is propagating 

vertically upward at infinite height, a value of p of co/a*, corresponds to a horizontal ray 

at infinite height, and a value of p of co/a« {[1 + az ] / [ 1 + az + AT/T.. ]} 1/2 corresponds to 

a ray with imaginary slope at infinity such that its turning point (point of horizontal slope) 
is located at the height z. Based on this description the waves group themselves into 
several different forms. 

The first group, 0 £ p £ p 0 = co/a„, {1 / [ 1 + AT/T.. ]} 1/2 , are waves with their turning 
points at or below the ground surface and thus are actually reflected at the surface. 

These waves plus their reflections constitute the first group. The wave with P = p 0 
grazes the surface and is the limiting ray between the reflected and refracted rays. The 
second group, p 0 < p < p 2 = co/a„, {[1 + az ] / [ 1 + az + AT/T.. ]} 1/2 , are waves with a 
turning point above the ground and below the observer height z. The waves in this 
group consist of those leaving the source in the range of angles described by (3.6) and 

their continuation after they have been refracted upward. The third group has p z £ p £ 

P s = co/a« {[1 + as ] / ( 1 + as + AT/T.. ]} 1/2 . These waves have their turning point above 
the observer and below the source. At the observers location these waves should yield 
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a exponentially decaying solution. The waves in this group consist of those leaving the 
source in the appropriate range of angles and their continuation following refraction 
upward. These three groups of rays can all be seen in a ray diagram for a point source 
and all initially are propagating downward from the source. In addition there are waves 

propagating upward initially. These are in the range 0 < p < P s but differ from the first 
groups in that they do not originate as downward propagating waves that are reflected 
or refracted upward. Thus the "reflection" coefficient is missing from these waves. 

In addition to the above groups that can be seen in a wave diagram for a point 
source, there are several types that are necessary for the superposition given by (3.3) 

where p ranges from zero to infinity, but do not occur in a ray diagram for a point 
source. Group five consists has P s < p < to/a^, these waves have the turning point 

above the source. In addition there are waves with co/a*, < p, these have no physical 
interpretation and correspond to complex angles at infinity. 

With these concepts let us proceed to the mathematical solution to the problem. 

The function g 3/2 given in (3.8) contains four branch points, two at P = ± P z and two at p 

i l 

= ± co/a*,. The negative branch points are not significant and will not be discussed. On 
the real p axis, for 0 < p < p z , g 3/2 (z, P) is real. For p 2 < p < co/a*, g 3/2 (z,p) is positive and 
imaginary, and for p > co/a.. g 3/2 (z,p) has a phase of -rc at p = co/a^and tends to a 
phase of -n/2 as p tends toward infinity. This is shown in Figure 3.1 . The branch line 
for g(z,p) = ( g 3/2 (z,p)) 2/3 can be chosen to be on the line where the phase of g 3/2 (z,p) is 
-k. This line extends from the first branch point at p z to the second at co/a*, and 
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encloses a small region above the positive real axis, see Figure 3.2. The branches 
chosen for g(z,P) yield a phase of zero for 0 < p < p 2 , tc/ 3 for p 2 < p < co/a«, and varying 

from -27C/3 to -rc/3 for p > co/a,.. 

Now by considering three cases the various forms of the solution can be 
obtained. These are shown in Figure 3.3. The first case is a wave with the turning 

point below the surface, 0 £ p £ p 0 . The second has the turning point below the source 

but above the surface, p 0 < p £ p s . In this case if the receiver is below the turning point 

then p > p 2 , if it is above the turning point then p < p 2 . The third case has p s <, p £ co/a« 
and the turning point is above the source. Again if the turning point is above the 

receiver then p > p 2 , if the receiver is above the turning point then p < P 2 . Note that 

these waves do not appear in a point source ray diagram but are needed to complete 
the solution. 

The solutions corresponding to the cases given above require the determination 
of the constants in (3.7). To do this a set of conditions are required. As result of the 
source terms in (3.4) the solutions separate into at least two forms, one for the region 
below the source and one above. The radiation condition, requiring outgoing waves in 
the limit as height tends to infinity, requires A to be zero above the turning point for 
z > s. At the source height, z = s, the solution must be continuous 

lim G(z,p) = lim G(z,p) 

Z— »S + Z— >8. 

(3.14) 

and must satisfy 
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. r 9G 9G. q 
lim [— ] - lim [—]=~ 

z-»s. dz z-»# + dZ 


(3.15) 

which is obtained by integrating (3.4) from z = s - e to z = s + e and taking the limit e-»0. 
At a turning point continuity is required, 


lim G(z,p) = lim G(z,p) 

Z— >Ztp . Z — >Zjp t 

(3.16) 

At the ground surface, z * 0, the required condition is the normal impedance condition 
(1 .4) which can now be expressed as 


G = 


i Z 9G 
(op dz 


(3.17) 

Using these conditions, equation (3.4), and the physical situations presented in 
Figure 3.2 the following solutions can be obtained. For z > s 


G = K h 2 (q(z,p)) [ fMnCs.p) + R 0 h 2 (ti(s,p)) ] 


for (o/a.. > p z > p s > p 0 > p > 0, 


(3.18) 
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G - K h 2 (T|(z,P)) [ Mti(s,P) + R, h 2 (ti(s,p)) ] 


(3.19) 


for oVa » > p 2 > p s > P > p 0 > 0, 


G = K h 2 (t|(z,P)) e in ^ R, [ MTtfs.p) + R 0 h 2 (r|(s,p)) ) 


for oVa » > P 2 > p > p s > p 0 > 0 and 


(3.20) 


G = K [ h 1 (ri(z,p)) e^+ R 0 h 2 (Ti(z,p)) ] 
Ri [ hi(Ti(s,p) + R 0 h 2 (n(s,P)) ] 


for > p > p 2 > p 5 > p 0 > 0. For z < s 


(3.21) 


G = K h 2 (Ti(s,P)) [ h^ruz.p) + Ro h 2 (rj(z,P)) ] 


for oVa^ > p 8 > p 2 > p 0 > p > 0, 


(3.22) 


G = K h 2 (n(s,p)) [ h^TKz.p) + R, h 2 (n(z,p)) ] 


for ca/a_ > p s > p 2 > p > p o > 0, 


(3.23) 
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G = K h 2 (r\(s,(J)) e** R, [ h^tUz.p) + R 0 h 2 (Ti(z,p)) ] 


for oVa.. > p s > p > p z > p 0 > o and 


(3.24) 


G = K [ Mrtfs.p)) e i7l/3 + R 0 h 2 (ii(s,p)) ] 
Ri[h 1 (TK2,p) + R 0 h 2 (n(z,P))] 

for oVa,, > p > P s > p z > P Q > 0. Where 


(3.25) 


K = K(g z (z,P), g z (s,p)) = 


12 i X 


f Vm z -P ) 1 >/ g z (s.P) 


(3.25) 


R o = ' 


2 (ri(0,P)) + i ¥ h 1 , (~n(0,P)) 
x h 2 (Ti(0,p)) + i v h 2 ’(ri(0 t p)) 


(3.26) 


x = a X - 


1 z gzzflf) 

2 P 0 a ~ g 2 (0,P) 


(3.27) 
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z 


and 


V = 


Po a „ 


(^) 3 g z (o,p) 


(3.28) 


.* 



(3.29) 

The solutions obtained above are for real values of p, but as discussed above p 
must be interpreted as a complex variable. The boundaries between solutions off the 

real axis must be chosen as the branch lines used for calculating the function g(z,P), 

g(s,p) and g(0,P) from g 3/2 (z,p), etc. as were discussed above. On crossing these 

branch lines it should be noted that the phase of g(z,p), g(s,P) and g(0,P) 

discontinuous^ jumps from -2 tc/ 3 to 2n / 3 and the phase of g z (z,P), g z (s,P) and g z (0,(3) 
increases by -2 tc/ 3 (since it contains the root of g in the denominator). 

Reference [8] presents some results for the modified Hankel functions that are 
useful for this type of behavior: 
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(3.30) 


and 

h.Oie 5 ) = -h 2 (n) 

(3.31) 

Using these relations the solution as given by equation (3.18) to (3.25) can be rewritten 
and the regions of validity determined. For z > s these are 


G = K(g z (z,P), g z (s,p)) {h 2 (ri(z,P)) [ h^TKs.0)) + Ro(g(0.p)) h 2 (n(s,p)) ]} 
in region A of p-space as given in Figure 3.4 


(3.32) 


G = K(g 2 (z,P), g 2 (s,P)) {h 2 (ri(z,p)) ( h 1 (ri(s,p)) + R o (g(0,p) 4**) h 2 (ri(s,P)) ]} 

(3.33) 


in region B 


G = K(g 2 (z,P)e ' 2n/3 , g 2 (s,p)) {h 2 (n(z,P)) [ h^p) e'™) 
+ R o (9(0.P) e> 2 ^ ) h 2 (n(s,P) e* 2^3 ) ]} 


(3.34) 


in region C 



G = K(g 2 (z,p)e* ‘ 2 **, g 2 (s,p) e* > 2** ) {h 2 (ri(z,p) e* *** ) [ h^nCs.p) e< ) 
+ R o (9(0.P) e' 2 * /3 ) h 2 (t|(s,P) e' 2x/3 ) ]} 


( 3 . 35 ) 


in region D. For z < s 


G = K(g 2 (z,p), g z (s,p)) h 2 (ri(s,p)) [ h 1 (ii(z,p)) + R 0 (g(0,p)) h 2 (Ti(z,p)) ) 


in region E 


( 3 . 36 ) 


G = K(g 2 (z,p), g 2 (s,p)) h 2 (Ti(s,p)) ( h^tKz.p)) + R„ (g(0,p) e* 2 ** ) h 2 (t|(z,p)) ] 


( 3 . 37 ) 


in region F 


G = K(g 2 (z,P), g 2 (s,P) e-^ ) h 2 (n(s,P)) [ h^rKz.p) e* 2 "*) 
+ Rq (9(0.P) e' 2 ** ) h 2 (r|(z,p) e' 2 **)] 


( 3 . 38 ) 


in region G and 


G = K(g 2 (z,p)e- * 2 */*, g 2 (s,p) e - i 2^ ) h 2 (Ti(s,p) e‘ 2*3 ) [ h^p) e* 2 */3 ) 
+ R 0 (g(O.p) e* 2n/ 3 ) h 2 (ri(z,p) e* 2Jt/3 ) ] 


( 3 . 39 ) 
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in region H. To obtain these results it is necessary to recognize that 


R,(g(0.W) - R o (g(0,p) e 3 ) 

(3.40) 

From these results and the description of the behavior of the function g(z,p) (and 

therefore r\) at the branch lines it should be clear that the solution is continuous at the 

branch lines even with g(z,p) being discontinuous. This transformed solution must 
now be inverted. 
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3.2 Lapse Results 


Two methods were used to approximately invert the Hankel transform contour 
integration which lead to the saddle point method in the insonified region of physical 
space and a FFT based numerical method. Neither of these methods will be described 
in detail as the basic method is well known in both cases and the specific application 
has been described in detail elsewhere. 

Both of these approaches are based on the concept that although the inversion 
integral (3.3) is defined as being carried out along the real axis, the residue theorem of 
complex variables [10] allows the path of integration to be changed provided that there 
are no poles of the integrand between the original and modified paths. If poles exist 
then additional terms must be included with the integral along the modified path. In the 
case of an isothermal atmosphere the additional term due to the pole leads to the 
surface wave term. In the lapse case the only possible pole is the due to the 
denominator of the term multiplying the upward going wave (the reflection coefficient in 
the case where a reflection occurs) being equal to zero. These case has not been 

completely examined but in the limit of AT equal to zero it reproduces the surface wave 
term. Thus one clearly expects to see a similar behavior in the case of weak lapse 
condition. This surface wave like-behavior has not been investigated beyond the point 
described above and has not been included in the results given below. 

3.2.1 Contour Integration-Saddle Point Method 

Integrals of the form 
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l(4)=Je' 5 Wl dp 
c 

(3.41) 

can be approximately evaluated by the saddle point method if the path of integration C 
is such that the ends of the path do not significantly contribute to the integral, £ is a 

large parameter, and f(p) has a point where the first derivative is zero. Complex 
variable theory indicates that a function can not have maximum in the region where it is 
analytic and a point where its first derivative is zero must be a saddle point [10]. At a 

saddle point if the path of integration follows the line of constant real part of f(p) then 
the imaginary part f(P) either increases or decreases at a maximum rate. If the path of 
integration follows the line of constant imaginary part of f(p) that passes through the 
saddle point then the real part of f(P) increases or decreases at a maximum rate. If we 

choose to follow a line of constant real part of f(p) through the saddle point in the 
direction such that the imaginary part increases at the maximum rate then the 
magnitude of the integrand decrease rapidly as we move away from the saddle point. If 

£, is large then the only significant part of the integral is near the saddle point and the 
integral can be approximated by using the first two non-zero terms in the Taylor series 

for f(p) yielding the well known results given in [1 1]. 

This method works well when a saddle point exists. However when one does not 
exist then an approximate integration can be carried out as described in detail by Ma 
[ 12 ]. 

To apply this method to the integrals given by (3.3) with (3.32) to (3.39) both the 
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Hankel functions and the modified Hankel functions must be replaced by their 
asymptotic expansions and the integral regrouped to extend from negative infinity to 
positive infinity. When this is done the result contains twenty terms since the integrand 

is different between in the various regions in (5-space and in each part of the G function 
contains two terms. Thus writing out only the terms of interest for z > s yields 



e 


3/2 

i{ P r + X[g (z,p) 


3/2 

g (s.P) l 


n/2} 


dP 


P. 

♦ J K 7 (z.P) 

K 


3/2 3/2 

e -'{Pr + Mg (z.P)-g (s,pn-*/2} 


dP 


+ ... 
P B 


~ 3/2 3/2 3/2 

fx ( 2 p)e* i{pr+M9 {z,p, + fl (0»]-«*) dp 


+ j K, 7 (z,p) e' 1 x 1 * o” *P> 1 “ * 1 dp 

Po 


+ ... 

(3.42) 
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/ 


and for z< s 


G(z.r) 


Po 

* 3/2 3/2 

+ Jc 6 (z,r)e -" |S '* lt9 M) - 9 

0 


3/2 ^2 

♦ Jc 7 (z.p)e " |l, * M9 ,s ' w ' 9 ftB1 " ,al dp 

Po 


Po 

- 3/2 3/2 3/2 

+ Jc i6 (z,|3)e l||l,<M9 (zW * 9 M) ' 29 «>J»]-*) d p 

0 

Pz ^ 

H. Jc i 7 (z.P) e ilPr + lt0 W>). 9 l W«*) d p 

Po 


+ ... 


( 3 . 43 ) 


To find the saddle points the arguement of of exponential term must be differentiated 
with respect to p and set equal to zero. On differentiating one finds 
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(*}) 

(3.44) 

where F(z,9(s)) is the function defined in (2.10) and used to describe the rays. 
Differentiating the arguments of the exponentials then yields equations (2.6) to (2.9), 
the equations defining the rays. Thus the saddle points associated with a particular 

point in physical space (in the insonified region) correspond to the values of p (or 0(s) 
where the two are related by (3.6)) that defines the two rays passing through the point 
in physical space. Just as the rays were interpreted as upward-going-direct waves, 
refracted waves, etc. the terms in (3.42) and (3.43) also have the same interpretations. 
Determination of the saddle point values then first requires determination of the types of 
waves present at a particular physical location and then solution of the appropriate two 
of equations (2.6) to (2.9). Once the location of the saddle point has been determined 
the classical results may be applied. A computer program for carrying out this 
procedure and the resulting equations to approximate the inversion integral have been 
given in detail by Cheng [13J and will not be repeated here. A typical result is shown in 
Figure 3.5. 

As the physical location of the receiver moves into the acoustic shadow real 

values of J5 or 0(s) cease to exist. Ma [12] has suggested an approximate approach 
which is also based on contour integration. In this approach the inversion integral 

between p 2 or p s and cc/a.. is carried out numerically and the remainder of the integral 

extending from negative infinity and to positive infinity are carried out in a manner 
similar to the saddle point method. The integral carried out numerically physically 
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represents the contributions of the exponentially decaying disturbances due to waves 
with turning points above the receiver. Ma [1 2] also incorporated into his program 
Cheng’s saddle point method for the insonified region with some changes. A typical 
result is shown in Figure 3.6. 

Both of these methods suffer from discontinuities as the receiver passes from a 
region where the receiver senses a direct and a reflected wave to one where a direct 
and a refracted wave would occur. This results not from the transformed solution which 
is continuous, but from the large argument approximation that must be made to obtain 
the saddle point form (3.41) and from the fact that the argument becomes zero at most 
of the boundaries. As a result of these intrinsic problems with the saddle point method 
a numerical approach was then applied. 

The saddle point method has the appeal of a physical interpretation of the 
mathematical steps and results. A purely numerical method loses that interpretation 
and the physical insight that comes from it. 

3.2.2 Numerical Integration Method 

The numerical approach used was developed by Richards and Attenborough [14] 
and was applied to the present case by Lloyd [1 5]. The method approximates the 
inversion integral by using a Fast Fourier Transform (FFT) algorithm. To obtain an 
integral suitable for the use of the FFT algorithm the Bessel function containing the 
horizontal distance dependence must be approximated by its asymptotic expansion. 
Three other modifications are then carried out. First, the integration path is modified to 
be above the real axis (Richards and Attenborough's original approach was to 
integrate below the axis but they also assumed e- kot .), this avoids the discontinuities at 
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the branch points but requires integration up the imaginary axis. Second, the integrand 
is modified to make the integral along the imaginary axis zero, but. as claimed by 
Richards and Attenborough, to not change the result. Finally an approximate term is 
added to account for the finite upper limit and to approximate the integral to infinity. 

This approach is described in detail by Lloyd [1 5] and Figure 3.7 is a typical 
result. Disadvantages of this approach are the large amount of computer time required 
and that a entire horizontal profile must be obtained at each receiver height. Thus to 
obtain a vertical profile many time consuming computer runs must be made and one 
point out of several thousand points is actually used from each run. This approach 
clearly does not contain the discontinuities present in the saddle point method. Figure 
3.8 compares the saddle point method and the purely numerical method. The 
agreement is excellent in the insonified region with the exception of the region very 
near the shadow boundary. The agreement is good in the initial sound level decrease 
as the shadow boundary is crossed but the saturation region deep in the shadow is not 
the same for the two methods. 

The numerical method often results in oscillations in the sound level at large 
distances from the source, this appears to be an artifact of the numerical inversion 
method and is dependent on the parameters of the inversion scheme. Also as very 
large distances are approached the calculated sound level often increases this is 
clearly due to the numerical inversion method. These points are further discussed by 
Lloyd [15] and a listing of Lloyd's program is given in Appendix A. 
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4.0 INVERSION SOLUTION 


The inversion solution, AT < 0, follows very closely along the lines of the lapse 
solution. The governing equation and boundary conditions, (1.2) to (1.4), are identical 
in the two cases and the approach using Hankel transforms is also the same. Again the 

solution requires analytic continuation off the real p axis. This has proved to be difficult 
and lead to several errors initially. (The solution given in the Sixth Semiannual Report 
[18] are incorrect.) 

Although the solution can be discussed in terms of rays and the closely related 
saddle point method, this approach has not been used to approximately evaluate the 
solution in the inversion case. In the lapse case only two rays, at most, pass through a 
given point. In the inversion case, at large distances from the source, many rays may 
pass through a given point due to the "trapping" effect of the inversion. Since the 
saddle point method requires all of these rays and their corresponding saddle points to 
be located, and this is the most difficult part of the method, the approach becomes 
impractical. Thus only the purely numerical method of inversion has been used. 

4.1 Transformation and approximate solution 

The time dependence is removed from (1.3) as in (3.1) and the resulting equation 
Hankel transformed using (3.2) to obtain (3.4). Again the location (in terms of the 

transform variable P) of the transition from oscillating to exponential behavior is given 

by (3.5). However, since AT/T.. < 0 the transition is at a value of p greater then co/a„. 
Using (3.6) and comparing results to those of Section 2.2 one can interpret the solution 
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in the range 0 < p < oo/a.. as representing the rays that will escape the trapping effect of 

the inversion. These rays either go directly from the source to infinite heights or go 
from the source to the ground where a reflection occurs and then go to infinite heights. 

In the range co/a.. < p < P s the solution represents rays that are trapped by the 

inversion. The transition given by (3.5) applies to rays in this range. Beyond this range 
the rays do not occur in a ray diagram. 

The solution of (3.4) can again be approximated by (3.7) and (3.9) through (3.1 1). 
Equation (3.8) must modified by a negative sign on the right hand side yielding 


£ 

9*(z.P) = -(1 + az +Y”) 



+ 



1 +Q 
1 -<D 


(4.1) 


This change is necessary since the region of oscillatory solution is below the turning 
point in the inversion case while it was above it in the lapse case (see Nayfeh [7]). In 

the inversion case AT/T.. < 0 and thus fl> > 1 for 0 < p < co/a.. It is convenient to note 
that we may rewrite the logarithm term in this case as 
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1 +p 

1 -o 


).Tanh'(i) + if 

<J> ^ 


(4.2) 


to clearly indicate the choice of the branch of logarithm, ln(-1 ) = +i n, the opposite of that 
in the lapse case. Thus, in the range 0 £ 0 < co/a^ g 3/2 (z,0) ranges from slightly below 
the negative real axis to infinity along the negative imaginary axis. For values of 0 
between oVa^ and 0 Z ranges along the positive real axis from infinity to zero. For real 0 

greater than 0 2 values of g 3/2 (z,0) are along the positive imaginary axis, see Figure 4.1. 
As in the lapse case the boundary between these regions are branch points of the 
function g 3/2 (z,0) with the branch lines extending downward from the branch points for 
Re(0) > 0 in 0-space. 

The argument of the Hankel functions involves g(z,0) = (g 3/2 (z,0)) 2/3 and again the 
branches must be chosen with care. For 0 < 0 < co/a« g(z,0) is chosen such that its 
phase ranges from slightly greater then zero (or 2k) at 0 = 0 to ic/3 as the branch point 

at 0 = io/a« is approached from values of 0 less then co/a,.. In this region the two 
modified Hankel functions have an oscillatory and exponential growth or decay 
behavior with with one (h,) representing upward traveling and decaying waves and 

the other (h 2 ) downward traveling, growing waves. In the range co/a« < 0 < 0 2 and 

g 3/2 (z,0) is real and positive. The function g(z,0) is chosen to be on the line with phase 
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-2x/3. In this case the modified Hankel function h 1 represents both upward and 
downward traveling waves, a standing wave-like phenomena, and h 2 a downward 
traveling wave. 

The branch line for g(z,p) = (g^z.P)) 2 / 3 needs to be defined to extend these 

solutions off of the real axis. Lines of constant phase of g^z.p) run froftn the first 
branch point to the second in a manner similar to the lapse case, but with the order of 
the branch points reversed. Figure 4.2 shows the behavior of these lines of constant 
phase. Again a line of constant phase is a convenient branch line. 

If the line where the phase of g 3/2 (z,P) equals -nJ2 is chosen as the branch line for 
g(z,P) = (g 3/2 (z.P)) 2/3 then the phase of g(z,p) can be made to agree with the desired 


values on the real axis as described above. In addition for real p and p 2 < p, g(z,p) 


has a constant phase of tc/ 3 with h^T^z.p)) having a decaying exponential behavior for 
increasing z and representing the contribution of waves with a turning point below the 
receiver's height to the total pressure field. 

Using the conditions given in (3.1 4) to (3.1 7) and the physical descriptions of the 
type of waves that occur in each situation the following solutions can be obtained. For 
z>s 

G = K [ h 1 (n(s,p)) + R 0 h 2 (q(s,p)) ] h 2 (ti(z,P)) 

(4.3) 


which is identical to (3. 1 8) for p 0 > P s > p z > co/a. > P > 0, 
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G = -K[ h 1 (ti(s,p)) / R 0 + h 2 (ri(s.P)) J h 1 (n(z.n 


(4.4) 


f0r P o > P» > P z >P>«ta..>0, 

.271 

h ■ 

G = - K e 3 [ h 1 (ii(s,p)) / R 0 + h 2 (Ti(s,p)) ] h 2 (T|(z,P) 


(4.5) 


for p o > p $ > p > P z > ca/a„, > 0 and 


G = K [ h/TKs.P)) + R 2 h 2 (ti(s,p)) ] h 2 (ri(z,P)) 


(4.6) 


(orp 0 > p > p 8 >P 2 > co/a^>0. Fors>z 


G = K [ h^Ttfz.p)) + R e h 2 (n(z.p)) ] h 2 (ri(s,p)) 


(4.7) 


the same as (3.22) for P 0 > P 2 > P s > a/a.. > p > 0, 


G = - K [ h 1 ('n(z.P)) / R 0 -*• h 2 (n(z.P)) ] h^nts.p) 


(4.8) 


forp o >P z > p 8 >P>w/a - >0, 
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(4.9) 


. 2 * 

G = - K e 3 | h,(r|(z,P)) / R o + h 2 (n(z.P)> ] hj(n(s.p) 


for J3 o > P 2 > p > p s >co/a 00 >Oand 


G = K [ h^(n(z,p)) + R 2 h 2 (n(z,p)) ] h 2 (fi(s,») 


(4.10) 


for P o > p > p 2 > p s > oVa.. >0. Here K and R 0 , and R, are defined by (3.25) and (3.26) 
and 


R 


2 



R 


o 



) 


(4.11) 


As discussed above these solutions must be continued off the real axis. As in the 
lapse case the boundaries are chosen as the branch lines for calculating g(z,p), g(s,p) 
and g(0,p) from g 3/2 (z,P), etc. On crossing these branch lines the phase of g(z,P), g(s.p) 

and g(0,p) jumps discontinuously from n/3 to - n and g 2 (z,p), g 2 (s,p) and g 2 (0,p) jumps 

discontinuously by 2rc/3. Using (3.30) and (3.31 ) solutions (4.3) through (4.10) can be 
continued off the axis as 
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G = K(g 2 (z,P), g z (s,p)) [ h 1 (ti(s,P)) + R o (g(0,P)) h 2 (Ti(s,p)) ] h 2 (r|(z,p)) 


in region A of Figure 4.3 


in region B, 


i 2* ,2ii . 4n 

G = K(g ; (z,P) e 3 , g 2 (s,P) e 3 ) [ h, (n(s,P) e 3 ) 


. An 


. An 


. An 


+ R o (g(0,P) e 3 ) h 2 (T!(s,p) e 3 )]hJii(z,P)e 3 ) 


(4.12) 


(4.13) 


- I* 


2 n 


.An 


G = K(g z (z.P). g 2 (s,p) e 3 ) { h,(Ti(s.p) e 3 ) 

. 4n .4 k 

+ R o (g(0.P) e 3 ) Mn(s.P) e 3 )] h,(t!(z,p)) 


in region C and 


.An 


(4.14) 


G = K(g z (z,P), g z (s,p)) [ h^rKs.p)) + R o (g(0,p) e ) h 2 (n(s,p» ] h 2 (ii(z,P)) 


(4.15) 


in region D. For z < s 


G = K(g z (z,p), g 2 (s,p» [ h/ntz.p)) + R e (g(o,p)) h 2 (n(z,p)) j h 2 (n(s,p)) 


(4.16) 
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in region E, 


I 2k . 2k . 4k 

G = K(g,(z,P) e 3 , g z (s,P) e 3 ) [ h,( n (z,P) e 3 ) 


. 4k .4k ,4k 

'T 'T 'T 

♦ R o (0(0.P) e ) h 2 to(z.P) e )] h 2 (Ti(s,P) e ) 


in region F, 


(4.17) 


.2 n . 4n 

G . K(g 2 (z,p) e 3 , g r (s,P) ) [ h, (n(z.p) e 3 ) 

. 4k 4k 

+ RJg(O.P) e 3 ) h,(n(s,P))l h.(r,(z,p) e 3 ) 


(4.18) 


in region G and 


. 4x 


G = K(g z (z,p), g z (®,p)) [ h 1 (ri(z,p)) + R o (g(0,p) e ) h 2 (ri(z.p)) ] h 2 (n(8,p)) 

(4.19) 


in region H. 

From these results and the description of the behavior of the function g(z,P) (and 

therefore t\) at the branch lines it should be clear that the solution is continuous at the 

branch lines even with g(z,p) being discontinuous. The transformed solution must now 
be inverted usirfg the numerical method developed by Richards and Attenborough [14]. 


* ( 
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4.2 Inversion Case Results 


As was discussed above the numerical method originally developed by Richards 
and Attenborough [14] was used to invert the Hankel transformed solution in the case 
of the inversion. This was due to the fact that many rays pass through a given point in 
physical space in the case of an inversion and the number of saddle points which exist 
equals the number of rays passing though that point. Since finding the saddle points is 
the most difficult and time consuming part of that method the approach appeared 
impractical in this case. The numerical method used is identical to that of the lapse 
case as described by Lloyd [1 5]. 

Only a limited number of cases have been run to date using the solution 
described in Section 4.1 and the numerical inversion technique. Figure 4.4 shows a 
typical case. The results generally show an interference pattern with 6 dB/doubling of 
distance decay out to distances of the order of thirty meters and a more complicated 
behavior beyond that distance but with no significant change in the rate of decay. This 
latter result is somewhat unexpected from qualitative arguments. Experimental data for 
propagation under inversion conditions is quite limited, with the data presented by 
Sutherland and Brown [16] being the major set. However, this set contains only seven 
measurements at a fixed height over a 675 meter distance. No direct comparisons 
have been made but the data also shows what appears to be a 6 dB/doubling of 
distance decay with some interference minima. Thus at least qualitatively the 
agreement appears good. 

Appendix B contains a listing of the program for the Inversion case. 
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5.0 CONCLUSIONS 


Approximate solutions of the Hankel transformed acoustic wave equation with a 
particular, realistic and well-developed vertical sound speed (or temperature) profile 
have been obtained for both the lapse and inversion cases. These solution are quite 
complex and exact inversion of the transformed solution does not appear possible. 

Both approximate inversion using contour integration and the saddle point approach 
and numerical inversion have been used to obtain the physical solution in the case of 
lapse conditions. Only numerical methods have been used with inversion condition. 

The lapse case shows the expected behavior: an interference pattern with a 6 
dB/doubling of distance decay within the shadow region; a rapid decrease in sound 
level in the vicinity of the geometric shadow boundary; and approximately a 6 
dB/doubling of distance decay well within the shadow region. Similar behaviors occur 
for both inversion methods but the contour integration - saddle point method yields and 
larger decrease in the sound level on passing into the shadow than the numerical 
integration technique. The origin of this difference has not been determined. The 
contour integration - saddle point method results appears to agree with the empirical 
model of Weiner and Keast [1 7] better then the results of the numerical inversion 
technique. Since the techniques are applied to the same approximate solution of the 
transformed acoustic wave equation the difference must result from the inversion 
techniques. The numerical technique also produces a weak interference-like behavior 
far into the shadow region. This appears to be artifact of the numerical method as is the 
increase in sound level that frequently occurs as the maximum distance for the 
inversion technique is approached. 

Agreement between the results and experimental data is fair within the shadow 
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boundary. The level is predicted well but the location of interference maxima and 
minima are not accurately predicted. This may be due to the poor fit of the temperature 
profile to the measured profile. No data appears to be available that both gives a 
temperature profiles and sound levels in the shadow region. 

The inversion case shows an interference pattern with a 6 dB/doubling of 
distance decay out to distances of the order of thirty meters for realistic temperature 
profiles. Beyond this distance the decay rate appears to remain nearly the same but 
the structure of minima and maxima becomes irregular. This tends to agree with a 
simple geometric argument since "trapped" rays start to reappear in a ray diagram at 
such distances. Little data is available for comparison in the inversion case. 
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7.0 UST QF SYMBOL S 
English 


a 

A 

B 

C 

° 1 . "^17 
D 
E 
f 



Sound speed. 

Function defined by (2.4) or constant in (3.7). 

Function defined by (2.5) or constant in (3.7). 

Function defined by (2.6). 

Constants. 

Function defined by (2.7). 

Function defined by (2.1 1). 

Arbitary function 
Function defined by (2.10). 

Function defined by (2.25). 

Function defined by (2.26). 

Function defined by (2.27). 

Function defined by (3.8). 

Hankel transform of G. 

Acoustic pressure with time dependence seperated out, see 
(3.1) 

Modified one-third order Hankel function of the first kind, see 

(312) : 

Modified one-third order Henkel function of the second kind, see 
(3.12). 

V- 1 

Intergal defined by (3.41). 

Zero order Bessel Function. 

Function defined by (3.25). 

Constants 

Acoustic pressure. 

Constant determining the strength of a point acoustic source. 
Horizontal distance from the source. 
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\ . 

I 

Horizontal distance from the source. 
Function defined by (3.26). 

Function defined by (3.29). 

Function defined by (4.11). 

Height of the point source above the ground. 
Time. 

Temperature. 

Height above the ground surface. 

Acoustic impedance of the ground surface. 


Scale factor for temperature, see (1 .2). 

Hankel transform variable replacing r, see (3.2). 
Function defined by (2.17). 

Delta function. 

Function defined by (3.10). 

Angle an acoustic ray make relative to horizontal. 
Limiting ray angle, see (2.28). 
a J(a) a). 

Arbitrary arguement 
Density of the air. 

Function defined by (3.27). 

Function defined by (2.12). 

Function defined by (3.9). 

Function defined by (3.29). 

Circular frequency. 
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Combinations 


AT Change in Temperature between the ground surface and 

far above the ground. 

Subscripts other than given above 

i 1 ,2 or 3. 

o Evaluated at the ground, or a reference value, 

tp Evaluated at a ray turning point, 

s Evaluated at the source height s. 

z Derivative with respect to height (g z or g Z2 ) or evaluated at the height z. 

«>«> Evaluated at infinite height. 
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Figure 1.1. Temperature as a function of height above the ground for different times of 
the day as determined by Best [1]. 


54 


302 


304 


306 306 

TE*>OUn« 00 


310 


312 


Figure 1 .2. The present model of temperature as a function of height and a set of 
observations by Butterworth {4]. 


55 




HEIGHT 



Figure 2.1 . The nomenclature used in defining the rays. 
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Figure 2.2. Acoustic rays for a lapse case with a = 1 .75 m' 1 and AT/T^ = 0.03 with a 
source at a height of 2 m. 
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Figure 2.3. Acoustic rays for an inversion case with a = 1 .75 m* 1 and AT/T^, ■ -0.03 
with a source at a height of 2 m. 
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Figure 3.1 . Path followed in complex g3/2 - space as the real part of p varies from zero 
to infinity and the imaginary part of p is constant for the lapse case. 


59 



Figure 3.2. Sketch of the location of the branch line used for calculating g (z,0) = 
( g 3/2 (z,p) j n the lapse case. 
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Figure 3.3. Sketch of the three types of physically occuring rays in the lapse case. 
Type 1 rays have their turning point below the ground surface. The turning point for 
type 2 rays is below the source and above the ground surface. Type 3 waves have a 
turning point above the source, this type of ray does not appear in a point source ray 
diagram but occurs in the superposition making up the inverse transform. 
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Figure 3.5. Typical results for a lapse case using the saddle point method only, from 
Cheng [13], as compared to the Weiner and Keast empirical model [1 7]. The solution 
extends only to the shadow boundary at about 68 meters. 
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Figure 3.6. Typical results for a lapse case using the combined saddle point-contour 
integration method. From Ma [12]. 
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Figure 3.7. Typical results for a lapse case using the numerical inversion technique. 
From Lloyd [15]. 
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Sound pressure level (dB) 



Figure 3.8. Comparison of the results of the saddle point-contour integration method 
and the numerical technique for a lapse case. 
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Figure 4.1. Path followed in complex g 3/2 - space as the real part of (3 varies from zero 
to infinity and the imaginary part of P is constant for the inversion case. 
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Figure 4.2. Sketch of the location of the branch line for calculating g (z,p) = 
( g3/2 (z,p) )2/3 j n the inversion case. 
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Figure 4.3. Sketch of the regions in complex P*space where the various forms of the 
solution are valid for the inversion case, the lines are branch lines for g (z,P), g (s,P) 
and g (0,P). 
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Figure 4.4. Typical results for an inversion case using the numerical inversion 
technique. 
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APPENDIX A * The Lapse Case Program 


The material below is from Lloyd [15] and both describes the program used to 
calculate the transformed solution and to carry out the inversion and presents a listing 
of that program. The program is named FFTPRESS and was written to be used on a 
VAX 750. Descriptions are given both of the subroutines comprising the program, of 
the input variables and a set of "helpful'hints" are given that may be useful in running 
the program. The intent of this section is not to describe in detail the operation of the 
program but to allow a somewhat experienced Fortran user to run the code as it was 
created. 

SUBROUTINES 

Input: Subroutine to input the necessary parameters to the main 
program. The following sentences summarize each of the input variables 
in the order they are requested. Tinf is the temperature at infinite height, 
normally 300 K. Tinf is used to calculate the speed of sound a. Dtot is the 
temperature change from infinity to the ground normalized by the 
temperature at infinity. Dtot is normally 0.025. Alpha is the term used in 
the temperature profile defining the altitude at which the temperature 
gradient becomes effective. Alpha is normally 2.5 (meters)" 1 . Splrefis 
the reference sound pressure level used in the calculations of sound 
pressure level in dB. Omega is the frequency of the sound source in 
rad/sec. Resistance is the flow resistance used in the Chessel model and 
is normally 300 cgs units. Zr is the height of the receiver in meters. Zs is 
the height of the source in meters. Alp is the amplitude of the imaginary 
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component in transformed space. Alp greater than zero is integration 
above the real axis. The product of Alp and the step size AK should not 
exceed 0.01 and may be much smaller. Alp is replace by (3 in the thesis. Ns 
is used in an analytical function that is subtracted from the sampled 
solution to null the effects of off axis integration. Ns is normally 3. If Ns 
is equal to zero then no subtraction occurs. Ns is replaced by 8 in the 
thesis. Me nonzero signals the inverse Hankel transofrm routine that the 
terms representing the integral extended to infinity are to be included in 
the inversion. Me is also the number of terms to be used and is normally 5. 
Me is replaced by M in the thesis. N1 is the number of points to be used. 

N1 depends on the maximum horizontal distance desired. N1 equal to 4096 
points is a common value. Np and N1 are used interchangeably. N1 must be 
equal to an integer power of two. Delbeta or delK are the step size in 
complex K space. In the program Delbeta and DelK are used 
interchangeably. Delbeta also depends on the maximum horizontal distance 
desired and also on the maximum Beta allowed. This maximum Beta is very 
near to omega divided by the speed of sound. Beta and K are used 
interchangeably in the program and thesis. 

Region: Subroutine used to determine which of 8 diffferent forms of 
the general solution are to be used. The selection depends on how the 
waves are interferring at that particular value of K. Region calls to 
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subroutine g32all to identify if the complex part of g 3/2 has changed sign. 
The sign change indicates a different set of rays are combining to form the 
solution. With each set of rays a different form of the solution is required. 

Hall: Subroutine to calculate the Hankel functions H-j and H2 and their 

derivatives in terms of the Airy functions, Al and Bl. Hall calls to Cgbair 
to get the Airy functions needed. 

Hall2: Subroutine to calculate only the Hankel functions. 

Cgbair: Subroutine to calculate the Airy functions. Cgbair uses either 
an asymptotic or a small argument approximation of Al and Bl depending on 
the value ofcomplex K. 

Gzalll: Subroutine to calculate the derivatives of the g 3/2 function. 
These values are used to compute the reflection coefficients. 

Gall: Calculates the g function needed to calculate the g 3/2 function. 

Dafb2: Subroutine modified from Attenborough and Richards to 
calculate the inverse Hankel transform using fast Fourier transforms. 

Dafb2 calls to subroutine Zeta to calculate the terms representing the 
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extension of the integral to infinity and to subroutine Fork used to perform 
the actual fast Fourier transform. 

Zeta: Subroutine to compute the value of the integral to infinity. 

Fork: An efficient fast Fourier transform taken from Attenborough 
and Richards program. 

VARIABLES 

Beta: Used interchangeably with K, both are the complex argument of 
the transformed solution. 

Tau: The term t used in the reflection coefficients developed by Van 
Moorhem. 

Sci: The term \y used in the reflection coefficient developed by Van 
Moorhem. 


Rq: The actual reflection coefficient. 


R-j : The modified reflection coefficient representing refraction. 
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En: The complex argument of the Hankel function. En is a function of 
K (or Beta) and the height of the source or receiver, whichever applies. 

Rlemda: The ratio of omega to speed times alpha. The reciprocal of 
the wave number. 

Zimped: The complex impedance of the ground normalized by the speed 
of sound and the density. 

Z2 and Z3: Heights of the receiver and source, respectively. 

Gbar (K) or Gbar(Beta): The sampled function to be inverted. 

Gbar(r): The inverted solution. The real space answer. 

G: The sound pressure level result. 

Rad: The horizontal distance of the present (Gbar (r) 


Rad2: The logarithm of Rad. 



HELPFUL NOTES 

1 . All units used here are examples only. The only requirement in 
operation is to use a consistent set of units. 

2. The output of sound pressure level and horizontal distance along 
with an echo of input parameters is written to file FOR044.dat in the 
present diretory. 

3. To plot the result at the University of Utah Mechanical Engineering 
Vax system type RUN PLOTPROMPT and the rest is interactive. 

4. Basic instructions for use on the University of Utah Mechanical 
Engineering Vax system are: 

a. Log on using normal sequency of user name and password. 

b. Type @Q to link all necessary files together. Instead of combining 
all files into a large file several small trackable files are used for ease of 
editing. 

c. Type RUN FFTPRESS to begin execution. 

d. Input the variables as requested by subroutine Input. 

e. At completion FFTPRESS will display FINALLY FINISHED. To plot the 
results type RUN PLOTPROMPT. This is a standard plotting program that 
uses the system subroutine Mgraph. The data file name is FOR044.DAT. 

The data file contains 2 columns. The fist column of the data is the 
logarithm of the horizontal distance. The second column of the data is the 
sound pressure level. 15 lines are used at the beginning of FOR044.DAT to 
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echo the input parameters, therefore tell PLOTPRMPT to begin accepting 
data at line sixteen of file FOR044.DAT. The plot will display on graphics 
terminals only. Mgraph asks if a hard copy is desired when crt plotting is 
finished. PLOTPROMT has autoscaling capability that can be turned on or 
off and offers many other self instructing options. Mgraph creates files 
named HPPLOT.HPL, however, it is recommended to change the name as soon 
as possible to avoid deletion of previous plots. If a hard copy plot is 
desired after exiting PLOTPROMT type PLOT then the file name, 
f. Log off with command LO. 


X 
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PROGRAM FFT PRESS 

MAIN ROUTINE TO DEVELOPE THE GBAR(b,Z) TO BE INVERTED * 
THE MAIN WILL INPUT THE NECESSARY CONSTANTS, CALCULATE BETA, * 
SELECT FORM OF SOLUTION CALL THE SUBROUTINES AND COMPUTE RO.TAU, * 
SCI, g, THE HANKEL FUNCTIONS AND FINALLY CALCULATE GBAR(b,Z). * 


SUBROUTINES: 

G32ALL: 

REGION: 

GZALL: 

GALL: 

HALL: 

HALL2: 

DAFB2: 


ZETA: 

FORK: 

INPUTS: 

TINF: 

DTOT: 

alfhaI 

SPEED: 

OMEGA: 

Z: 

S: 

SPIREP: 

RESISTANCE: 

ALP: 

NS: 

ME: 

N1: 

DELEETA: 


FINDS THE FUNCTION g‘3/2(b,z) USED TO DETERMINE 
WHICH REGICN OF SPACE PRESENT BETA IS IN. 

GIVEN g*3/2 DETERMINES WHICH REGION OF 
SPACE PRESENT BETA IN. 

FINDS gz(*,Z)=dg/dZ AND THE OTHER DERIVATIVES 
USED TO FIND K,TAU,SCI. AISO RETURNS ETA(b,Z) 
THAT IS USD AS THE ARGUMENT FOR HANKEL 
FUNCTIONS. GZALL CALLS TO GALL TO GET g. 

CALCULATES -g FUNCTION. 

CALCULATES THE HANKEL FUNCTIONS FROM AIRY 
FUNCTIONS. CALLS CGBAIR TO GET AIRY FUNCTIONS 

CALCULATES ONLY THE HANKEL FUNCTIONS NOT THE 
DERIVATIVES AS HALL DOES. CALLS CGBAIR ALSO. 

GIVEN GBAR(b,Z) USES METHOD DEVELOPED BY 
RICHARDS AND ATTENBOROUGH TO PERFORM THE 
HANKEL INVERSION. USES SUBROUTINES ZETA, FORK. 
DAFB2 MAKES SEVERAL CORRECTIONS TO A GENERAL 
FFT. THE STANDARD CODE IS TAKEN FROM RICHARDS 
AND ATTENBOROUGH PROGRAM. 

FORMS THE SUM OF ME TERMS WHICH APPROXIMATES 
GBAR(b,Z) TO INFINITY IN THE BETA SPACE. 

A VERY FAST FFT USED TO PERFORM ACTUAL 
INVERSION OF GBAR(b,Z) FROM THE BETA SPACE. 

THE TEMPERATURE AT VERY LARGE Z 

THE DELTA_T/T PARAMETER REPRESENTING THE 
TEMPERATURE GRADIENT. 

PARAMETER USED IN DEFINITION OF TEMPERATURE 
GRADIENT. 

SPEED OF SOUND AT T(INF). 

FREQUENCY OF SOUND IN RADIANS /SBC. 

FIXED DISTANCE TO THE OBSERVER Z2 IN PROGRAM. 

FIXED DISTANCE TO THE SOURCE Z3 IN PROGRAM. 

REFERENCE SOUND PRESSURE LEVEL USED TO COMPUTE 

THE SOURCE STRENGH Q. 

GROUND RESISTANCE IN THE CHESSELL MODEL. 

THE TERM USED TO INTEGRATE OFF REAL BETA AXIS. 

PARAMETER IN THE ANALYTICAL FUNCTION IN DAFB2. 

PARAMETER TO PRODUCE SUM TO INFINITY IN DAFB2. 

SIZE OF ARRAY TO BE INVERTED. 

STEP SIZE USED FOR BETA ALSO DELK IN DAFB2. 


NOTE: K AND BETA AND DELK AND DELBETA ARE USED INTERCHANGEABLY 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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VARIABLES: 

BETA: INDEPENDENT VARIABLE IN TRANSFORMED SPACE. 

TAU: TAU DEFINED IN PAPER BY VAN MOORHEN. 

SCI: SCI DEFINED IN PAPER BY VAN MOORHEN. 

RO: RO CONSTANT DEFINED IN VAN MOORHEN. 

R1 : R1 CONSTANT DEFINED IN VAN MOORHEN. 

EN: HANKEL FUNCTION ARGUMENT. 

RLENDA : OMEGA/ ( ALPHA *SPEED ) . 

ERO: BRANCH CUTS IN THE BETA SPACE ASSOCIATED WITH 

BRZ: THE SQUARE ROOT AND 3/2 POWER FUNCTIONS 

BRS: IN g AND g THREE_HAiF. 

BRW: BRANCH CUT AT OMEGA/SPEED. 

Z IMPED: GROUND IMPEDENCE NORMALIZED BY DENSITY AND SPEED. 

GZ: dg/dZ FROM SZALL1 . 

GZZ: d2g/dZ2 FROM GZALL1 

OTHER DERIVATIVES f>ER THIS NOTATION 
Z1 : REFERENCE DISTANCE O.O. 

Z2: Z AS ABOVE. 

Z3: S AS ABOVE. 

MU IK K K tHHUHHHHUHHU M SHHHHHHHHHHHHHHHt m H UHK OKK I tiKKKIlKKIlKKKKKKlIIlKlIKKK 

IMPLICIT DOUBLE PRECISI0N(A-H,0-Z) 

INTEGER IRBGION 
COMMON /lNTEG/'NS,MMl 
COMMON /AFB2IN/ ALP, DELBETA .. 

COMMON /C0NSTANT1 / SPEED, OMEGA 
COMMON /C0NSTANT2/ ALPHA, DTOT 
COMMON /C0NSTANT3/RLEMDA , Q 
COMMON /C0NSTANT4/CMPI , CK , RO , R1 
COMMON /HEIGHT/ Z1 ,Z2,Z3 
COMMON /CETA/Z1EN,Z2EN,Z3EN 
COMMON /BRANCH/BRO , BRS , BRZ , BRW 
COMPLEX* 1 6 BETA, GZ, GZZ 
COMPLEX* 1 6 H2,H21 
COMPLEX* 1 6 EN,Z1EN,Z2EN,Z3EN 
C0MPIiK*l6 HI ,H1 1 
COMPLEX*! 6 CK, TAU, SC I, Z IMPED 
COMPLEX* 16 DUM1 , DUM2.R0, R1 ,CMPI 
C0MPLEX*16 GEAR (32768) 

CALL INPUT (TINF.SPIREF, RESISTANCE) 

Q= . 00002*4 .*3.141 5926*4 . 67DO* ( 1 0 . ** ( SPLREF/20. DO) ) 

SPEED=DSQRT(1 .4DO*287.DO*TINF) 

PRINT *, ’THE FOLLOWING IS AN ECHO OF THE INPUT ’ 

PRINT *, ’IN THE FOLLOWING ORDER ALP DELTA ME NP DELBETA’ 

PRINT *, ’SPEED OMEGA ALPHA DTOT Z1 Z2 Z3 RESISTANCE’ 

PRINT *, ’TINF.SPIREF’ 

PRINT *, ’ ’ ! SKIP A LINE 

PRINT *, ’THESE VALUES ARE ALSO WRITTEN TO FILE 44 ’ 

PRINT *, ALP, NS, ME, N1 , DELBETA, SPEED, OMEGA, ALPHA, DTOT 
PRINT *,Z1 ,Z2,Z3, RESISTANCE, TINF.SPIREF 

WRITE(44,*) ’ECHO OF INPUT ALP NS ME DELBETA SPEED OMEGA ALPHA 
& DTOT Z1 Z2 Z3 RESISTANCE TINF SPIREF’ 
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VRITE(44,*) ALP, NS, ME, HI , DELBETA , TINF, SPIi^EP, OMEGA , ALPHA , 

& DT0T,Z1 ,Z2,Z3, RESISTANCE, Q, SPEED 
RLEMDA=OMEGA/ ( SPEED*ALPHA ) 

PI =3. 141 5926 
CMPI=(0.0,1.D0) 

C FUNCTIONS TO COMPUTE THE COMPLEX IMPEDENCE FROM CHESSELL MODEL 
FR=0MEGA/2. DO/PI 

RATIO=FR/RESISTANCE , V 

R=1 .+9-08D0*RATI0**(-.75D0) 

X=>-11 .9DO*RATIO**(-.73DO) 

ZIMPED=DCMPLX(R,X) 

BRO=OMBOA/SPEED*DSQRI £ 1 .DO/(1 .DO V ALPHA*Z1+DTOT)) 
BRZ=OMBGA/SPEED*DSQRT , ( ( 1 .DOfALPHA*Z2)/(1 . DO+ALPHA*Z2+DTOT ) ) 
BRS=OMBGA/SPEED*DSQRT ( £,1 . DOf ALPHA*Z3 ) /( 1 . DO+ALPHA*Z>DTOT ) ) 
BRW=OMBGA/ SPEED 

PRINT *, • ' * 

PRINT *, ' ' r 

PRINT *, 'THE BRANCH CUTS ARE' ,BRO,BRZ,BRS,BRW 
WRITE(44,*) 'THE BRANCH CUTS ARE’ ,BRO,BRZ,BRS,BRW 
DO 1 ,1=1 ,N1 

BETA =DCMPLX ( DFLOAT ( 1-1 ), (-ALP))*( DELBETA) 

IP (ABS(BETA) .LB. .OOOOOOOl) THEN 
BETA=DCMPLX [• OOOOOOT DO, (-ALP) ) *DELBETA 
END IF 

CALL GZALL1(Z1,BETA,GZ,GZZ,JN) 

Z1EN=EN 
DUM1 =GZ 
DUM2=GZZ 4 

CALL GZALL1 (Z2,BETA,GZ,GZZ,FN) 

Z2IN=EN 

CK=Q/1 2.D0/CMPI/(RLEMDA**(2.D0/3.D0) )*1 .0/(CDSQRT(GZ) ) 

CALL GZALL1 (Z3,BETA,GZ,GZZ,EN) 

Z3EN=EN 

CK=CK*1 .DO/(CDSQRT(GZ)) 

TAU=ALPHA*RLEMDA-CMPI/2 . D0*ZIMPED*DUM2/DUM1 

SCI =4 IMPED* ( ( 3 • /2 . ) ** ( 2 . DO/3 • DO ) ) * ( ( RLEMDA ) « ( 2 . DO/3 . DO ) ) *DUM 1 

CALL HALL(Z1EN,H2,H21 ,H1 ,H11 ) 

DUM1 =TAU*H1 4CMPI*SCI*H1 1 
DUM2=TAU*H2+CMPI*SCI*H21 
R0=-DUM1 /DUM2 
DUM1 =CMPI*PI/6.D0 

R1=(CDEXP(-DUM1 )*(CMPI*R0))/((CDEXP(DUM1 )**2.D0)+(R0**2.D0) ) 
CALL REGI ON ( BETA , IREGION ) 

GO T0(10, 20, 30, 40, 50, 60, 70, 80), IREGION 
C ********** REGION 1 BEGINS HERE FOR Z2>Z3 OR Z>S ********** 

10 CONTINUE 

CALL HALL2(Z3EN,H2,H1 ) 

DUM1 =H1 -fR0*H2 

CALL HALL2(Z2IN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION 2 BEGINS HERE FOR Z2>Z3 OR Z>S ********** 
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20 CONTINUE 

CALL HALL2(Z3EN,H2,H1 ) 

DUM1 =H1 +R1 *H2 

CALL HALL2(Z2EN,H2,H1 ) 

GEAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ****** ** ** REGION 3 BEGINS HERE FOR Z2>Z3 OR Z>S ** * ** ***** 
30 CONTINUE 

CALL HALL2(Z3EN,H2 f H1 ) 

DUM1 = (HI +R0*H2)*R1 *CDEXP(CMPI*Pl/3. DO) 

CALL HALL2(Z2EN,H2,H1 ) 

GEAR ( I ) =CK*H2*DUM1 
GOTO 500 • * 

C ■» » « » »»» REGION 4* BEGINS HERE FOR Z2>Z3 OR Z>S ^ « «««» 
40 CONTINUE 4 \ 

CALL HALL2(Z3EN,H2,H1 ) 

DUM1 =CDEXP ( CMPI*Pl/3 .DO ) *R1 *(H1 +R0*H2) 

CALL HALL2(Z2EW,H2,H1 ) 

GBAR(I)=CK*DUM1*(H2+(CDEXP(CMPI*PI/3.D0)*H1)) 

GOTO 500 

C ********** region 5 BEGINS HERE FOR Z2<Z3 OR Z<S ********** 
50 CONTINUE 

CALL HALL2(Z2EN,H2,Ht ) 

DUM1=H1+R0*HT 

CALL HALL2(Z3EN,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION € BEGINS HERE FOR Z2<Z3 OR Z<S ********** 
60 CONTINUE 

CALL HALL2(Z2EN,H2,H1 ) 

DUM1 =H1 +R1 *H2 

CALL HALL2(Z3M,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION 1 BEGINS HERE FOR Z2<Z3 OR Z<S ********** 
70 CONTINUE 

DUM1 =CDEXP(CMPI*PI/3.D0)*R1 
CALL HALL2(Z2EN,H2,H1 ) 

DUM1 =DUM1 *(H1 +R0*H2) 

CALL HALL2(Z3*N,H2,H1 ) 

GBAR ( I ) =CK*H2*DUM1 
GOTO 500 

C ********** REGION 8 BEGINS HERE FOR Z2<Z3 OR Z<S ********** 
80 CONTINUE 

DUM1 =CDEXP ( CMPI *PI /3 . DO ) *R1 
CALL HALL2(Z2EN,H2,H1) 

DUM1=DUM1*(H1+R0*H2) 

CALL HALL2(Z3EN,H2,H1 ) 

GBAR(I ) =CK*(H1 *CDEXP(CMPI*PI/3. D0)+H2)*DUM1 
GOTO 500 

C END OF GBAR (BETA) CALCULATIONS BASED ON REGIONS DETERMINED 
C MULTIPLY BY BETA ONLY TO MATCH VAN MOORHEM DEFINITION TO 
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C RICHARDS AND ATTENBOROUGH. 

500 CONTINUE 

GBAR(I)=GBAR(I)*BETA 

1 CONTINUE 

PRINT ' LAST BETA EXCUTED IS', BETA 
CALL DAFB2(GBAR) 

DO 2,1=1 ,Nl/2 

RAD=2 . DO*PI *DFLOAT ( 1-1 )/(DFL0AT(N1 )*DELHETA) 
IP(RAD.LE.0.0 ) THEN 
GOTO 5 
END IP 

RAD2=DL0G1 0 ( RAD ) 

G=20. DO*DLOG1 0( ABS (£BAR( I ) ) ) 5 - 

WRITE (44, 9093) RAD2,G 
5 CONTINUE 

2 CONTINUE 

PRINT *, 'FINALLY FINISHED' 

9093 F0RMAT(5X,3G18.8) f 

STOP * 

END 


SUBROUTINE INPUT (TINF,SPIREF*RESISTANCE) 


C THE PURPOSE OF THIS ROUTINE IS TO INPUT ALL NE CESSA RY PARAMETERS TO * 
C THE MAIN ROUTINE. THE DEFINITION OF EACH PARAMEOTt WILL BE DEFINED * 
C AT ITS RESPECTIVE INPUT. * 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON /INTEG/ NS,ME,N1 
COMMON /AFB2IN/ ALP,DELBETA 
COMMON /CONSTANT 1 / SPEED, OMEGA 
COMMON /C0NSTANT2/ ALPHA, DTOT 
COMMON /C0NSTANT3/ RLEMDA,Q 
COMMON /C0NSTANT4/ CMPI,CK,R0,R1 
COMMON /HEIGHT/ El ,Z2,Z3 
COMMON /CETA/ Z1 EN,Z2EN,Z3EN 
COMMON /BRANCH/ BRO, BRS, BRZ , BRW 
WRITE (6,899) 

WRITE(6,900) 

READ *, TINF 
WRITE (6,901) 

READ *, DTOT 


WRITE (6,902) 
READ *, ALPHA 
WRITE (6,903) 
READ *, SPIREF 
WRITE (6,904) 
READ *, OMEGA 
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WRITE (6,905) 

READ *, RESISTANCE 
WRITE (6,906) 

READ *, Z2,Z3 
Z1=0.D0 
WRITE (6,899) 

WRITE (6,907) 

READ *, ALP 

ALP=-ALP ! NECESSARY DUE TO PREVIOUS DEFINITION USED 
WRITE (6,908) 

READ *, NS 
WRITE (6,909) 

READ *, ME v 

WRITE (6,910) . 

READ *, N1 
WRITE (6,91 1 ) 

READ \ DELBETA . . 

899 PORMAT( 'O', 20X.1H ,///////////////) ! CLEARS THE SCREEN 

900 FORMAT ('O’ ,///,lOX,3$J THE FOLLOWING ARE PHYSICAL VARIABLES. 

& ,//,12X,45H INPUT THE TEMPERATURE AT Z EQUAL INFINITY. 

& ,/) 

901 FORMAT (X,12X,45H INPUT THE TEMPERATURE CHANGE DELTA T/T(INF). 

* /) 

902 FORMAT (X,12X,46H INPUT THE TEMPERATURE PROFILE CONSTANT ALPHA. 

& ,/) 

905 FORMAT (X,12X,46H INPUT THE REFERENCE SOUND PRESSURE LEVEL. 

& ,/) 

904 FORMAT (X,12X,46H INPUT THE ANGULAR FREQUENCY IN RADIANS /SEC. 

& ,/) * 

905 FORMAT (X,12X,47H INPUT THE CHESSELL MODEL FLOW RESISTANCE (cgs) 

& ,/) 

906 FORMAT (X, 12X.46H INPUT RECEIVER AND SOURCE HEIGHTS SEPARATED 

& ,/,25H BY A COMMA.,/) 

907 FORMAT (///,10X,37H THE FOLLOWING ARE NUMERIC VARIABLES. 

& ,//,12X,39H INPUT ALP. THE IMAGINARY PART OP BETA. 

& , ,/,12X,45H NOTE POSITIVE VALUES ARE ABOVE THE AXIS. ,/) 

906 FORMAT (X,T2X,46H INPUT THE DELTA PARAMETER USED IN INTEGRATION 
& ,/) 

909 FORMAT (X,12X,46H INPUT THE NUMBER OF TERMS IN SUM TO INFINITY. 

& ,/) 

910 FORMAT (X,12X,46H INPUT THE NUMBER OF POINTS TO BE USED. 

& ,/) 

911 FORMAT (X, 1 2X.46H INPUT THE STEP SIZE IN BETA (OR K) SPACE. 

& ,/) 

RETURN 

END 


SUBROUTINE REGION ( BETA, IREGION) 
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C SUBROUTINE TO DETERMINE WHICH OP 8 REGIONS * 

C THE EVALUATION IS TO TAKE PUCE IN. * 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

INTEGER DOUBLE PRECISION K1 
COMPLEX*1 6 BETA.G32C1 ,G32C2,G32C3 
COMPLEX*! 6 G32PI,G32PI2,G32PI3 
COMMON /INTEG/NS,ME,N1 
COMMON /AFB2IN/ALP , DELBETA 
COMMON /CONSTANT! /SPEED, OMEGA 
COMMON /CONSTANT2/ ALPHA , DTOT 
COMMON / CONSTANT3/RLEMDA , Q v 

COMMON /CONST ANT4/CMPI ,CK,RO,R1 
COMMON /HEIGHT/ Z1 ,Z2,Z3 
COMMON /CETA/Z1EN,Z2EN,Z3EN 
COMMON /BRANCH/ERO, BRS r BRZ , BRW * 

INTEGER IREGION 
CALL G32ALL(Z1 , fcETA,G32C1 ) 

CALL G32ALL( Z2 , BETA , G32C2 ) 

CALL G32ALL(Z3,BETA,G32C3) 

C IN THE EVENT THAT BRZ OR BRS AND BRW ARE VERY CLOSE THE FOLLOWING 
C IS USED TO INSURE PROPER , REGION IS CHOSEN. W1 AND W2 ARE SIMPLE 
C WEIGHT FACTORS TCL DETERMINE WHICH DIRECTION TO EVALUATE REGIONS. 

C BY WEIGHTING THE SELECTION DEPENDING HOW CLOSE BETA IS TO BRW 
C THE REGIONS ARE FOUND IN THE REVERSE ORDER. 

P0SMAX=MAX(Z2,Z3) 

Wl=1.0 
V2=1 .0 

DUD=OMEGA/SPEED 

DUD=DUD*W1 +DUD*W2*DSQRT ( ( 1 . Of ALPHA*POSMAX ) / ( 1 +ALPHA*POSMAX+DTOT ) ) 
DUD=DUD/2 
IF(Z2.LT.Z3) THEN 
GOTO 120 
END IF 

C FOR **(Z2>Z3 OR Z>S)** THE FOLLOWING IDENTIFY REGIONS OF SPACE 
A11=DIMAG(G32C1 ) 

AI2=DIMAG(G32C3) 

AI 3=DIMAG ( G32C2 ) 

IF ( REAL ( BETA ) .GT. DUD) THEN 
GOTO 110 
END IF 

IF(AI3.GT.0.0) THEN 
IREGION =4 
GOTO 150 
END IF 

IF(AI2.GT.0.0) THEN 
IRE)GI0N=3 
GOTO 150 
END IF 

IF(AI1 .GT.O.O) THEN 
IREGI0N=2 
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GOTO 150 
END IP 
IRBGI0N=1 
GOTO 150 
110 CONTINUE 

IF( AI1.LT. 0.0) THEN 
IREGIQN=1 
GOTO 150 
END IP 

IP(AI2.IT.0.0) THEN 
IREGIQN=2 
GOTO 150 

END IF 1 

IF(AI3.IT.0.0> THEN 
IREGIQN=3 
GOTO 150 * 

END IF . t . 

IREGI0N=4 

GOTO 150 ■/ . 

C FQR**(Z2<Z3 OR Z<S)** THE FOLLOWING IDENTIFY REGIONS OF SPACE 
1 20 CONTINUE 

AI 1 =DIMAG ( G32C 1 ) 

AI 2=DIMAG ( G32C2 ) 

AI 3=DIMAG ( G32C3 ) . 

IF ( REAL ( BETA ).GT. DUD) THEN 
GOTO 130 
END IF 

IF (AI3.GT.0.(T) THEN 
IREGI0N=8 •• • 

GOTO 150 
END IF 

IF (AI2.GT.0.0) THEN 
IREGI0N=7 
GOTO 150 
END O’ 

IF "'(All .gt:o.o) then 
IRBGI0N=6 
GOTO 150 
END IF 
IREGI0N=5 
GOTO 150 
130 CONTINUE 

IF (AI1.LT. 0.0) THEN 
IREGI0N=5 
GOTO 150 
BJD IP 

IF (AI2.LT. 0.0) TH0J 
IREGI0N=6 
GOTO 150 
END IP 

IF (AI3-LT. 0.0) THEN 
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IREGION-7 
GOTO 150 
0JD IP 
IREGIQN-8 
GOTO 150 
150 CONTINUE 
RETURN 
DO) 

C 

C 

C 


C • * 

SUBROUTINE G32 ALL ( Z , BETA , G32C ) 



IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*16 BETA, PHI, SQRT1 ,SQRT2,F00,G32C,M0DL0G 

C0MPLEX*16 SI ,A,B,S 

COMMON /C0NSTANT1 /SPEED, OMEGA 

COMMON / C0NSTANT2/ALPHA , DTOT 

COMMON /BRANCH/BRO, BRS , BRZ , BRV 

AA=1 . +ALPHA*Z+DTOT 

B=1 . -( ( SPEED*BETA /OMEGA ) * ( SPEED*BETA/OMBGA ) ) 

SI =SQRT1 (BETA,Z) 

S=SQRT2(BETA) 

PHI -SI /(S*DSQRT(AA) ) 

IF(DT0T.EQ..0.0R.CDABS(1 . -PHI ).LT. ID-8) THEN 
FOO-.O 

ELSE 

POO=MODLOG( ( 1 .+PHI )/( 1 .-PHI ) ) 

ENDIF 

G32C=#SQR? CAA )*S1 5*DT0T*P00/S 

RETURN 
END 
C 
C 

C BEGIN OP FUNCTIONS USED ABOVE. 

C 

C SQRT1 AND SQRT2 ARE FUNCTIONS TO CALCULATE SORT (BETA) GIVEN * 

C DESIRED BRANCH CUTS AND DIRECTION -(-IMAGINARY OR -IMAGINARY. * 

FUNCTION SQRT1 (BETA.Z) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX* 16 BETA.AA.SQRTI 
COMMON /CONSTANTI /SPEED, OMEGA 
COMMON /CON STANT2 /ALPHA , DTOT 

AA=1 .-( ( (SPEED/OMEGA)*BETA)*( (SPEED/OMBGA)*BETA) ) 
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EB=1 .+ALPHA*Z+DTOT 

BRAN1 =( OMBG A/S PEED )*DSQRT ( ( 1 .+ALPHA*Z )/( 1 . +ALPHA*Z+DTOT ) ) 
IP ( (DREAL(BETA) .GE. BRAN1 ) .AND. 

1 (DIMAG (BETA) .LE. O.O)) THEN 
SQRT1 =-CDSQRT ( AA*BB-DTOT ) 

ELSE 

SQRT1 =CDSQRT ( AA*BB-DTOT ) 

end if 

RETURN 

END 


C 

c 

FUNCTION SQRT2(BETA) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX* 1 6 BETA,AX,SQRT2 
COMMON /ERANCH/BRO , BRS , ERZ , BRW 
COMMON /C0NSTANT1 /SPEED, OMEGA 
COMMON / C0NSTANT2 /ALPHA , DTOT 
AA=1 .-( ( ( SPEED/OMEGA )*BETA)*( (SPEED/OMBGA)*BETA) ) 
IF ((DREAL(BETA) .GE. BRW) .AND. 

1 (DIMAG (BETA) .I£. .0).) THEN 
SQRT2=-CDSQRT ( AA ) 

ELSE 

SQRT2=CDSQRT(AA) 

END IF 
RETURN 
END 


C 

c 

c 

c 

C FUNCTION MQpLOG COMPUTES THE LOG OF BETA GIVEN DIRECTION AND * 

C LOCATION OF BRANCH CUTS. * 

FUNCTION MODLOG(QUAN) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX*! 6 QUAN,MODLOG 

IF ( (DREAL(QUAN) .LE. 0.0) .AND. ( DIMAG (QU AN) 

1 .GE. 0.0)) THEN 

MODLOG=CDLOG ( QUAN )+DCMPLX (0. 0, -2*3 .141 5927 ) 

MODLOG=C DLOG ( QUAN ) 

END IF 
RETURN 
END 


C 

C 

C 

C 


END OF FUNCTIONS USED ABOVE. 
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SUBROUTINE HALL(Z,H2,H21 ,H1 ,H1 1 ) 

HALL USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER 
HANKEL FUNCTIONS FROM AIRY FUNCTIONS. 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX*1 6 Z,AI ,BI, AIP,BIP,K,KS,H1 ,H2,H1 1 ,H21 
C0MPLEX*16 ARG.CI 
CI= DCMPLX(O.DO,1 .DO) 

PI= 3.H1592654DO 

ARG= DCMPLX(O.DO,-PI/6.DO) 'v 

K= (l2.DO)**(1.DO/6.DCr)*CDEXP(ARG) 

KS= DCONJG(K) * 

CALL CGBAIR (-Z.AI.BI^AIPiBIP) 

H1= K*(AI-CI*Bl) 

H2= KS*(Al4CI*BI); > 

H1 1* -K*( AIP-CI *BIP) 

H21 = -KS* ( AIP+CI*BIP ) * 

RETURN 

END 


SUBROUTINE CGBAI R(Z,AI,BI,A IP rBIP) 

C CALCULATE AIRY FUNCTIONS FOR COMPLEX*! 6 ARGUMENT * 

C REP. HANDBOOK OF MATHEMATICAL JUNCTIONS, AERAMOWITZ AND STBGUN. * 
C ENTRY: • 

C CALCULATE ARGUMENT(Z) AND ABSOLUTE VALUE(Z) * 

C IF /Z / LT 6 * 

C THEN USE EQS. 10.4.2 THRU 10.4*5 FOR AI,BI,AIP,BIP * 

C 10 ELSE IF ARG(Z) LT PI/3 * 

C THEN CALCULATE ZETA(Z) * 

C USE EQSf 10.4.59, 10.4.61, 10.4.63, 10.4.66 FOR AI,BI,AIP,BIP* 
C 20 ELSE CALCULATE ZETA(-Z) * 

C USE EQS. 10.4.60, 10.4.62, 10.4.64, 10.4.67 FOR AI,BI,AIP,BIP* 

C ENDIF • 

C ENDIF * 

C EXIT * 

C END * 

Q m**«t*******HH****<***»»t****W*«)HHHh»*«*«««*»H(***«fHHH(«t 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 Z,A1,BI,AIP,BIP,ZETA,CZETA,Z14,SUM1 ,SUM2,SUM3,SUM4, 
1 ZETAP,FACT1 ,FACT2,SN,CS,FTERM,FPTERM,GTERM,GPTERM,F,FP # G,GP,Z3 
C0MPLEX*16 VZETA,VZETAP 
DIMENSION C(21),D(21) 

DATA Cl ,C2, PIRT, PI4/.355O20O539DO, . 25881 94038D0, 1 .772453851D0, 
+ . 7853981 635DO/ 

DATA C/1 .DO, .069 4444 4444 4 4 44D0, 

+ .0371 35487654321 DO, .0379930591 27800DO, 
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1 .0576491 9041 2669D0,. 11 609906402551 DO, 

+ . 291 591 59923074D0, . 8776669695099800, 

2 3.0794530501 731 DO, 1 2.341 573332345DO, 

+ 55- 62278536591 4D0, 278. 465O0O77759DO, 

3 1533.1 6943201 27D0, 9207. 2065997258D0, 

+ 59892. 51 356587500,41 9524.8751 1653DO, 

4 31 48257. 41 78666D0, 251 9891 9-871 601 DO, 

+ 21 4288036 • 96366D0, 1 929375549 ♦ 1 823DQ, 

5 18335766937.889D0/ 

DATA D/1 .DO, 

+ -.097222222222221 DO, -.0438850308641 97D0,-.042462830789894D0, 

1 -.0626621 63492031 DO, -. 1 241 0589602727D0, -. 306253764901 OTDO, 

2 -.92047999241 291 DO, - 3 .<21 04935846485D0, -1 2. 807 29306073 5D0, 

3 -57. 50830351391 IDO, -287.0332371 0920D0,-1 576.3573033370D0, 

4 -9446 . 3548230953D0, -6^335 • 70666384700,-428952 . 40040004DO, 

5 -321 4536 . 521 4006D0, -25697908. 383909D0, -21 8293420.8321 4D0, 

6 -1 963523788. 9909D0, r 1 8643931 088. 105D0/ 

ABSZ=ABS(Z) 

IP(ABSZ.BQ.O) GO TO 3. 

IF(ABS(DIMAG(Z)).LE.1.D-12.AND.DREAL(Z).I/T.0.D0) GO TO 5 
ARGZ=ATAN2(DIMAG(Z ) , DREAL(Z) ) 

GO TO 4 

3 ARGZ=O.DO 
GO TO 4 

5 ARGZ=3 • 1 41 5926535898D0 

4 CONTINUE 

IF(ABSZ.GT.6.D0) GO TO 10 
C ASCENDING SERIES 

C EQS. 10.4.2,10.4.3 

15 CONTINUE 
Z3=Z**3 
FTiRM=1 .DO 
FPTERM=Z*Z/2 . DO 
GTERM=Z 
GPTERM=1 .DO 
GLIM=1 4 D-1 3*ABSZ 
F=FTERW 
FP-FPTERK 
G=GTERM 
GP=GPTERM 

KKKT=100 ! ADJUST KKKT TO INSURE CONVERGENCE IP NECESSARY 
DO 1 1=1 ,KKKT 
13=3*1 

FTERM=FTERM*Z3/ ( (13-1 .D0)*I3) 

FPTERM=FPTERM*Z3/ ( 13* ( 13+2 . DO ) ) 

GTERM=GTERM*Z3/(I3*(I3+1 .do)) 

GPTERM=GPTERM*Z3/( (I3-2.D0)*I3) 

F=F+FTERM 

FP=FP+FPTERM 

G=G+GTERM 

GP=GP+GPTERM 

IP(ABS(GTERM) .LE.GLIM) GO TO 2 
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1 CONTINUE 
PRINT 6000, Z 

6000 PQRMAT(/' Z='2E14-5, ' ERROR IN CGBAIR, NONCONVERGENCE* ) 

2 AI=C1*F-C2*G 
AIP=C1 *FP-C2«GP 

BI=1 .732O5O806DO*(C1*F4O2«G) 

BIP=1 .752050808D0*(C1 *FP+C2*GP) 

GO TO 9999 

C ASYMPTOTIC EXPANSIONS FOR /Z/ LARGE 

10 SIGN=1 .DO 
SUM1=0.D0 

SUM2=0.D0 . 

SUM3=0.D0 * 

SUM4=0-D0 * 

PIBY3=3 • 1 41 5926DO/3 . DO -• 

IF(ABS(ARGZ) .GE.PIBY3) GO TO 20 , 

/ARG(Z)/ LE PI/3 • t 

EOS. 10.4.59, 10.4.61, 10.4.63, 10.4.66 

ZETA=CZETA ( ABSZ , ARGZ ) * 

DO 11 1=1 ,12 
K=I-1 

ZETAP=ZETA**K 
SUM1 =SUM1 +SIGN*C ( I ) /ZETAP 
SUM2=SUM2+SIGN*D ( I ) /ZETAP 

SUM3=SUM3+C(I)/ZETAP 

SUM4=SUM4+D ( I ) /ZETAP 

1 1 SIGN=-SIGN ... , 

Z1 4=ABSZ** . 25DO*DCMPLX(COS (ARGZ/4. DO) , SIN(ARGZ/4 .DO) ) 
FACT1 = . 5D0»EXP ( -ZETA ) / ( PIRT*Z1 4 ) 

FACT2= . 5D0*EXP (-ZETA ) *Z1 4/PIRT 
AI=FACT1*SUM1 
AIP=-FACT2*SUM2 
FACT1 =EXP ( ZETA ) / ( PIRT*Z 1 4 ) 

FACT2=EXP(ZETA)*Z1 4/PIRT 
BI=FACT1 *SUM3 
BIP=FACT2*SUM4 
GO TO 9999 

/ARG(Z)/ GT PI/3 NOTE CHANGE ABOVE 
BQS. 10-4-60, 10.4.62, 10.4.64, 10.4.67 

20 CONTINUE 

ARGZ=ATAN2(-DIMAG(Z) , -DREAL(Z) ) 

ZETA=CZETA ( ABSZ , ARGZ ) 

VZETA=1 .DO/ZETA 
LLL=10 

DO 21 1=1, LLL 
K2=(I-1 )*2 
J=K2+1 

VZETAP=VZETA**K2 
SUM1 =SUM1 +SIGN*C( J )*VZETAP 
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SUM2=SUM2+(SIGN*C( J+1 )*VZETAP*VZETA) 
SUM3=SUM>(SIGN*D(J)*VZETAP) 

SUM4=SUM4+ ( SIGN*D (J+1 )*VZETAP*VZETA) 

21 SIGN=-SIGN 

Z1 4=ABSZ** . 25DO*DCMPLX ( COS (ARGZ/4 . DO) , SIN ( ARGZ/4 . DO) ) 
FACT1 =1 .D0/(PIRT*Z14) 

FACT2=Z1 4/PIRT 
SN=SIN ( ZETA+PI 4 ) 

CS =C0S ( ZETA+PI 4 ) 

AI=FACT1 *(SN*SUM1 -CS*SUM2) 

AIP=-FACT2*( CS*SUM3+SN*SUM4 ) 

BI =FACT 1 * ( CS*SUM1 +SN*SUM2 ) V 

BIP=FACT2*( SN*SUM>CS*SUM4 ) 

9999 RETURN 

IND '* 

« 

.v * • 

BEGIN OF FUNCTIONS USEE ABOVE 


FUNCTION CZETA ( ABSZ , ARGZ ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMPLEX*! 6 CZETA 
ARG=ARGZ*1 . 5DO " 

CZETA=(ABSZ**1 . 5D0)*DCMPLX(C05(ARG) ,SIN(ARG) )*.66666666666667DO 

RETURN 

IND 


0JD OF FUNCTIONS USED ABOVE. 


SUBROUTINE GZALL1 ( Z , BETA , GZ , GZZ , IN ) 

C GZALL1 CALCULATES ALL THE PARTIAL DERIVATIVES * 

C OF THE g fUNCTIQN. * 

IMPLICIT DOUBLE PRBCISION(A-H,0-Z) 

COMPLEX *16 BETA,GZ,GZZ,G,GB,GEB,SQRT1 ,EN 
COMMON /CDNSTANT1 /SPEED, OMEGA 
COMMON /C0NSTANT2/ ALPHA , DTOT 
CALL GALL(Z,BETA,G) 

A=1 . +ALPHA*Z+DTOT 
BRAN=OMEGA/SPEED*SQRT ( ( A-DTOT ) /A ) 

IF (DREAL(G).LE. ,O.AND.DIMAG(G).GE.O. ) THEN 
SI=-1 . 

TgT.cre 

SI=1. 

h3^ fDJP 

GZ=SI*2.*ALPHA*SQRT1 (BETA,Z)/(3.*CDSQRT(G*A) ) 

C=2 . *ALPHA**3 • DO*DTOT/( 9*A**2 . DO ) 

GZZ=C/(GZ*G)-. 5*GZ**2. DO/G 
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T=(3. *0MBGA/(2.*ALPHA*SPEED) )**(2.D0/3.D0) 
EU=T*G 
RETURN 
END 
C 
C 

c 

c 

SUBROUTINE GALL (Z, BET A, G) 

C 

C GALL EVALUATES THE g FUNCTION * 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX *16 G,BETA,G32 
COMMON /CONSTANT 1 /SPEED-; OMEGA 
COMMON / C0NSTAMF2 /ALPHA , DTOT 
CALL G32ALL(Z,BEIA*G32)- 
G=CDEXP( 2 . /3 . *CDL0G(G32 ) ) 

RETURN ; 

END 

C 

C 

C ' . 

C - 

SUBROUTINE HALL2(Z,H2,H1 ) 

C HALL2 USES SUBROUTINE CGBAIR TO CALCULATE 1/3 ORDER * 

C HANKEL FUNCTIONS FROM AIRY FUNCTIONS. NOT THE * 

C DERIVATIVES AS HALE DOES. * 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMPLEX*1 6 Z,AI,BI,AIP,BIP,K,KS,H1 ,H2,H1 1 ,H21 
COMPLEX* 1 6 ARG,CI 
C0MPLEX*16 BETA 
CI= DCMPLX(O.DO,1 .DO) 

PI= 3.L41592654DO 

ARG= "DCMPLX ( 0 . DO, -PI /6. DO) 

K= (12.DO)**(1 .DO/6.DO)*CDEXP(ARG) 

KS= DCONJG(K) 

CALL CGBAIR(-Z,AI,BI,AIP,BIP) 

HI = K*(AI-CI*BI) 

H2= KS*(AI+CI*BI) 

RETURN 

END 

C 

C 

C 

C 

SUBROUTINE DAFB2(F) 

C SUBROUTINE TO ACCURATELY DO THE HANKEL TRANSFORM OF THE SOUND * 
C PRESSURE LEVEL. * 
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C F(NP)=GBAR(NP) MUST BE SAMPLED AT NP POINTS WITH K=(N-1, ALP) * 

C ALP REPRESENTS THE DISTANCE ABOVE THE REAL AXIS THE FUNCTION WILL • 

C BE INTECRATED. • 

C NS IS A PARAMETER REPRESENTING ADDITION OP AN ANALYTICAL FUNCTION * 

C TO F(NP) • 

C M IS THE NUMBER OP TERMS USED TO APPROXIMATE P(NP) TO INFINITY * 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /INTEG/NS,ME,N1 
COMMON / AFB2 IN /ALP , DELBETA 
COMPLEX*l6 F(N1 ),CF,CARG,SUM,FNP,CMPI,D1 
NP=N1 \ 

CMPI=DCMPLX(O.DO, 1 .DO) 

deuc=delbeta 

PI=3 - 1 41 5926DO . '• 

C SUBTRACT THE ANALYTICAL FUNCTION IP- NS > ZERO 
C ADJUSTING THE SUBTRACTION MULTIPLIER CP 
IP(NS.LE.O) GOTO »1 1 
CF=DCMPLX(O.DO,O.DO) ; 

IP (ALP.E3Q.O.O) THEN 
CP=DPLOAT (NP) /DFLOAT (NS ) 

CF=CF*F(2) 

END IF 

IP (ALP. NE. O.O) THEN 

CF=CMPI *DPLOAT ( NP ) *F ( 1 )/(DFK>AT(NS)*ALP) 

END IF 

C SUBTRACT THE ANALYTICAL FUNCTION IF NS>0 
DO 10,1=1 ,NP i 
D1 =DCMPLX( DFLOAT (1-1 ) , (-ALP) ) 

CARG=DFL0AT(NS)*(-D1 ) /DFLOAT (NP) 

P(I)=F(I)-CP*(1 .DO-CDEXP(CARG) ) 

10 CONTINUE 

1 1 CONTINUE 

IF(ALP.EQ.O.O) F(1 )=DCMPLX(O.DO,O.DO) 

FNP=F(NP) 

DO 12fI=2,NP 

D1 =DCMPLX ( DFLOAT (1-1 ) , (-ALP) ) 

F(I)=F(I)/(CDSQRT(D1 )) 

12 CONTINUE 

IP(ALP.NE.O.O) F(1 )=F(l )/CDSQRT((-CMPI)*(ALP)) 

C ADD TERMS TO INFINITY IP ME>0 
IF(ME.LT.1 ) GOTO 20 
DO 15,1=1 ,NP 

D1 =DCMPLX( DFLOAT ( 1-1 ) , ( -ALP ) ) 

CF=D1 /DFLOAT (NP) 

CALL ZETA(NP,ME,CF,SUM) 

F(I)=F(I)+FNP # SUM 
15 CONTINUE 

20 CONTINUE 

C DO THE FFT 

CALL FORK(NP,F, 1 ) 

C ADD ALTERNATE TERMS TO GIVE NP/2 SAMPLES TRANSFORMED 
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CF=DEUC*DFLOAT ( NP )* ( CDSQRT ( -CMPI ) ) / ( 2 . DO*PI ) 
DO 25,I=2,NP/2 

A=DEXP ( DFLOAT ( 1-1 )*(ALP)*2.DO*PI/DPLOAT(NP)) 
F(I)=A # F(I )+CMPI*F(NP-I+2)/A 
F(I)=F(I)*CF/DSQRT(DFL0AT(I-1 )) 

25 CONTINUE 
RETURN 
PHD 


C 

C 

C 

C *■ v 

SUBROUTINE FORK^LX.CX.SIGNI) . 

C HUM 

C A FAST FFT GIVEN BY J.F. CLAMBOUT, "FUNDAMENTALS OF GEOPHYSICAL * 
C DATA PROCESSING" . , . * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 CX(LX),CARG,CW,CTEMP,CI,DUM1,CMPI 
INTEGER SIGNI 


10 

20 


30 


J=1 

CMPI=DCMPLX(0.D0, 1 .DO) 

PI =3- 141 5926 „ 

SC=DSQRT ( 1 . DO/DFLOAT ( LX ) ) 
DO 50, 1=1, LX 
IF(I.GT.J) GOTO 10 
CTMP=CX(J)*SC 
CX(J)=CX(I)*SC 
CX(I)=CTEMP 
M=LX/2 

IF(J.LE.M) GOTO 30 

J=J-M 

M=M/2 

IF(M.GE.I) GOTO 20 
. * J=J+M 
L=1 - 4 - 


40 ISTEP=2*L 

DO 50,M=1 ,L 

CARG=CMPI *PI *DBLE ( SIGNI )*DBLE( (M-1 ) )/DBLE(L) 
CW=CDEXP(CARG) 

DO 50,I=M,LX,ISTEP 
CTEMP=CV*CX(I+L) 

CX(I+L)=CX(I)-CTEMP 
50 CX(I )=CX(I)+CTEMP 

L=ISTEP 

IF(L.LT.LX) GOTO 40 

RETURN 

END 
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SUBROUTINE ZETA(NP,M,A,SUM) 

f W w w W n * W 

SUBROUTINE TO ADD THE NECESSARY TERMS TO 
EXPRESS INVERT I BLE FUNCTION TO INFINITY. 

WILL USE DOUBLE PRECISION. 

SUM=SUM OF 1/(NP\5)*1/((J+A)\5) FOR J=1 TO 
INFINITY MINUS SOME CONSTANT WHICH IS 
INDEPENDENT OF A. 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

C0MPLEX*16 A,SUM,D2 
D2=DCMPLX(DFLOAT(M),O.DO) , 

SUM=1 .DO/(M+A) . , v 

SUM=2.DO*(DSQRT(DFLOAT<M))-1 .DO/CDSQRT(SUM)) 

$ -0.5 # CDSQRT(SUM)*(1.O+SUM*(1. 0/1 2.0-SUM*SUM/l 92.0)) 

DO 10, J=1 ,M 

SUM=SUM+1 . DO/CDSQST ( J+A ) 

10 CONTINUE 

SUM=SUM/DSQRT ( DFBOAT (NP ) ) 

RETURN 

END 


APPENDIX B - The Inversion Case Program 


The code for the lapse case follows. This code is extremely similar to the lapse 
case and the subroutine names and functions, variable names and hints are identical 
or at least very similar to those in the lapse case. 
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PROGRAM MAIM 


ORIGINAL PAGE IS 

OF POOR QUALITY 




c 

c 

c 

c 


c 

c 


c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 


Halo prograa: will call the subroutine* tot 

1) Input the program parameters. 

2) Fine tha points in tha g*3/Z w ,fu«ctlon there 
region changes occur. 

3) 3uxlo tha a* t r i i of values obtained by torching 
slightly above tha real axis# 

4) p er f ore tha Hankal Transform on the tatrlx 
obtained in step three. 

5) Print the results. 


implicit g: u?le 

COMMON / y 4lN/ 
CC*MCh /CONST/ 
CCMMCN /C0NST1/ 
COMMON /C0NST2/ 
common / C ON S T 3 / 

CCMMCN /REGION/ 
COMMON /ST AT:/ 
COMMON / I N T £ 0 / 
COMMON /iFfLIN/ 


PRECISION (A-N,0-Z> 
f 

Cl, PI 

SPSEC, OMEGA 

ALPMA/OTCT 

Z/ZREF/S 

R$1,RS2/RZ1,RZ2,R01,R02 
PCS, RESISTANCE 
M/Mg/NOPTS 
hGMT, CELESTA 


C*.MPl£X*1o r ( 327ed),CI,GlESS,ROOT,G32 

c: = :cmplx< o.ccei.oo) 

= 0ATANi(C.C0#-1.3C> 
p A I N T * / * 3 : s ' , p 2 


CALL INPUT<0tL5cTA,HGhT,RT1,RT2) 
A N 0 = -PI/2.00 


PRINT *,'06TEkP:n:*G region CHANGE coordinates* 
call REGION. FIN0(ANG/RT1/RT2) 


PRINT *,'8i:L0tNG MATRIX' 

CALL SUIL0M4TRIX (CSL3ETA, MGH T ) 


PRINT ** 'OC INC HANK cL TRANSFORM* 
CALL HANK5L 


print *, # Prirting results to F0RC45.0AT* 
CALL PR INTA Li ( OEtBETA, NOPTS ) 


PRINT COMPlETc •** ^ 7 —f ' ~ 

STOP ' 
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r* r. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ORIGINAL PAGE IS 

OF POOR QUALITY 



SUBROUTINE INPUT (0EL8 ETA/ MGHT/RT1/*T2> 


this is the input section op the program 
ill parameters necessary por the operation 
CP the PROGRAM hill be input in this routine 


c 

c 

IMPLICIT OCOELS PRECISION (A-H/O-I) 

CCMMCN /CONST/ CI/PI 

/ccnstt/ speeC/Omega 

CCMMCN /CONSTE/ JLPHA/CTCT 
SvHMCK /CON >T 3 / I/IRSF/S 
COMMON /STATE/ :CE/RcSISTANCE 
COMMON /inte:/ * S / m£ / NOP T s 
CCMPLEX‘1o Cl 

c 

•RITE (c / i OC ) 

PRINT »/'In s LT s^TSICAL PARAMETERS: Cataetor haight/ Sourca haight/' 
PRIsT •/ ' *nd Angular Fraquancy' 

READ 

PRINT 'INPLT STATE P a R am E T E R S : T aao gradiant/Taeparatura profile' 
print »/* constant and Tapp, at infinity' 

REAC • / 0 T OT / A l P h A / T IN F 

PRINT ‘/'Input tH# cFassall rodtl r LCH RESISTANCE* 

R * AC ‘/RESISTANCE 

1 PRINT ‘/'Input tra NUMBER OF OATA PCINTS' 

R E AC ‘/NCPTS 

TST ‘CL OG COFLOAT ( NO PTS>) /CLOG (2. COT 
IrCCAESCTST-CNlNT (TST)T .GT. 1.0-8) THEN 
NCPTS*ONINT<E.CO“CNINTCTST)J 

PRINT ‘/'Invalid input Caust fca a poaar of 2) TRT AGAIN* 

P* ENT • /' (a*. '/NCPTS/')^ 

GOTO 1 
END IF 

PRINT */'IN?LT nANKEL TRANSFORM PARAMETERS: NS/ ME and Intagratlon 

i naight' 

R E AC */NS/M£/ mGMT 
JREF * C.CO 

SPEEC * OSCRT (ACT .SCO « TINP) 

RT1 * OMEGA/SPEEC 

RTE ‘ RT1.0SCRTC1. 00/C1. CC+OTOT)) 

CEL3ETA « RTE • T.01CC /NCPTS 

NRITE <6/300) " 

8GC F0RNATC‘C'/2CX/1M//////////////////////W/> f elaars thp acraan 
PRINT ‘/'Rafaranca hght *•*/» t* 

PRINT ‘/'Oatactor haightf»*VZ 
PRINT ‘/'Sourca Haight *'/S 
PRINT ‘/'Fracuancy Cr ad)* - */ OMEGA 
PRINT ‘/'spaad of sound s*'/S*ff# “ "• 

PRINT •/ * Taap Gradient »*/OTOT V 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


PRIM *, # ProfiIa (alpna) « # ,UPMA 

PRINT *, # T infinity » # ,Tinf 

PRINT *, # Cbassal rasist.* # ,R6SISTANCC 
D RINT * 0 # Nua • points * # ,NCPT$ 

PRINT *,'0al data * # *0El%ET4 

PRINT *,'- # 

PRINT *,' Karkal' 

PRINT NS **/NS 

PRINT*,' * E 5',P£ 

print */ ' 4 Ic * # / h G MT 

PRINT prat ion h;nt=',HGHT*OEL!STA 

PRINT *, 
n c T L R N 

; N 0 


Sw5RCwTIN= at IL CP -TRIX(JCI!ET4, UP) 


Tms sL-rcutina is jsad to dataraina tha function 
*V m * e n locrnti's notis. It sill fca callad by 
tna Tam prp^rem ana will return tba aatri* able* 
<al« ce irvartap fcy tha bankal progna. 


IMPLICIT : CJ5LE 
common /main/ 

Cw'^C;, /CONST/ 

common / const w 

: CM HCN /CONST! / 
CO V NON /CONST 3/ 
COMMON /RSGICN/ 
CjM^CN /STATE/ 
common / int eg / 


PR-CISIONCA — ,C-2) 

GEAR 
Cl, PI 

crtEC, OMEGA 
ALPHA, 3TCT 
i/ifiEF/S 

cS1,*?S2,fi21,*Z2,R01,ff02 
SCE, RESISTANCE 
KS#«E#NOPTS 


C0MPlcX*1o S2AR<:i76E)y3ETA 


CCMPIEX*16 GZ_I,GZ,S,GZ_OGZZ_OGZZ,GZ 
C0MPlEX*16 iTA.Z^ETA. S/ETA. 0 


0 C M P L E X * 1 o rt_N:*r2_NZ*Hl_NS/H2_NS 
C0MPL£X*16 rtl.NO H2.NC/ nll^NOJ H21.HO 
C0M?LEX*16 CCNST,R,E/CI/T1,T2/T3 
COMPLEX*! 0 TAU#PSI#ZIMP 
CCMPLSX*16 SPGICsSRGZSsSRGZsSRGS 

integer iregion 

E = CDEXPC- PI * Cl / 3.00) 

rR = CMEGA/ C2.0OPI) 

PATIC = C R/R£$ISTANCE 


00 


X 

Z I HP 
Q 
Q 

RLHOA 

RLHOA23 


1 .CC*9 .CS OC* RATIO* *(-.7500) 
-11.90C*RATI0**(-.73DQ> 

OCPPLX ( 00#X ) 

1.C0 

O.COOOZ*4.CO*PI*4.67DO*(1 0 .00**C SPIRE P/ 20.00 ) ) 
CHEGA/ (SPEEO*ALPMA) 

R L * C A •* (2.00/3.00 
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•• PROGRAM BEGINS HERE ** 


C 

C 

c 

c 

c 

c 


c L 

c 


hL 


origin; 

OF POOR QUALITY 


co ic,z»i,n:pts 


BETA * OCPPLX (OFLCAT(I-I)/ (ALP))*(DELBETA) 

CALL GZ*LLl (ZREF,dETA,GZ.C/GZZ.O,ETA_0> 

CALL GZALL1 (Z/3=TA,GZ.Z,GZZ,ETA_Z> 

CALL GZALLl (S/CcTA/GZ.S^GZZ/ETA^S) 

CALL HALL (STA_C*h2.NC/H21 _NC/Hl _NC#Ml 1 _NC# IREGO) 

CALL ntLcZ(£TA - S/n2.K-/n1.NS/IRrGS) 

CALL *ALL2 (ETA.Z,H2.NZ,-i1.NZ,IREGZ> 

REGION (£ETA,IS£GICM 

:al = i.:c - <s?c=:/c??ga ♦ ci /2 * zihp • gzz.o/gz.O) 

3 s: * : : k c * c : * s ? - £ : / c v = c a 

'Si = *s: * <(3. • rlpoa>**<2.do/3.oo)> * gz.c 

s^s: s ::scat(c:.z> 

S R 3 3 s E 3 S C * T ( G Z _ S ) 

::>.gt = ; • eeta / (’z*ci*(rlmca23>*srgz*srgs> 

C(ir^(rt2.Nu)) .EC. C.OC .ANC. (A8S(F21.NQ> . EQ • 0.00) ThEN 
• s cc m plxci .cie^c.co 

GOTO cl 

i\Z I r 

: r 7A'j • h2_.SC ♦ PSI * F21.N0 

; * -CTAU * - 1 _ \C ♦ =»SI * FII^NO) / R 

SC T Z i 1 » l * 1 * 4 » £ / t / 7 / $ ) * I R c G I C N 


c 



I 

3 

I 

» 

4 

I 

i 

5 


i 

t 



* RESIGN 1 * 

Ga^C) = C C N 5 T * (M.N$ ♦ >-2_N$*R) * N2.NZ 

SCTO 10 

* REGION 2 * 

GdAR(I) « CC NS T * (hi .NZ ♦ H2.NZ*R> * H2.NS 
GOTO 10 

* REGION 3 * 

GdARCZ) *• (CONST /(E*c)> • (H1_NSZR ♦ H2.NS ) * H2_NZ 
GOTO 10 

* REGION 4 * 

GSARO =-CCNST,* (h1_N$/R ♦ H2. N S ) * Hi. HI 
GOTO 10 

* REGION 5 * 

G 3 A R ( I ) *-CCNST * CM.NZ/R ♦ M2_NZ) * Hi. NS 
GOTO 10 

** REGION 6 ** 

G3AR(I) » CONST * (HI. NS ♦ <<R-E)/R)*E * H2.NS > a H2.M2 
GOTO 10 

* REGION 7 * ‘ 

GdAR(I) » CONST • (hi .NZ ♦ <<R-EJ/*)*E * H2.N2) * H2.NS 
GOTO 10 


• REGION 6 * 


100 


or# 0000000 000000 000 


ORIGINAL PAGE IS 

OF POOR QUALITY 


G B AR ( I ) « CCNST « e • (hl.NZ/R ♦ M2.N2) * M2.NS 


CONTINUE 

RETURN 

END 


Subroutine gjjali C/ e eta,g J2C) 


S32AUL CALCULATES $3/2 FUNCTION 


» .«c3ifi«d to tF# or ancF ooint of tFa LOG function 

• is a • i 0 * * " 1 # cositiv* p» «1 axis. F«b 20/ 86 

» •noa.ifna #i*in 1C/12/BE for mocifitd proflraa 


implicit ocueli precision <a— «/C-2> 

CCM?LEX*l6 SETA/P N I / S C S T1 / SC B T 2/ FOO/G32C/NOOLOG 

complexes si, A/S/S/C: 

CCMMCN /CONST/ CI/PI 
CCMMCN /CCNST1 / SPEEC/CN50A 
CC-MCN /CONST!/ AL?rA/CTCT 
A A * 1 .00 ♦ ALP*-A»i ♦ OTCT 

i * 1.;0 - ( ( SX£E2*=tTA/0HEGA) * <SPEEC*6£TA/0NEGA> ) 
SI * S;«T1 (EETA/C) 
s * sc-tkeetaj 

?nl x SI /(S*CSC«T(AA)> 

I-CGTCT .EC. ,C . CR . CO AB S ( 1 . -onl ) .UT. 1.0-8) THEN 

'03 s 0.0 

a «.S E 

f00« .5 * CCLCGC <1.fPHl)/(1.-PNl) ) 

: NO I F 

032C ■ - OS C fi T (AA) ♦ SI ♦ CTOT * FOO/S 

R c TURN 

ENO 


6cGIN OF FUNCTIONS L SEC ABCVE. 



SCRT 1 A NO SCAT 2 ARE FjK'CTICNS TC CALCULATE SCRTIBETA) GIVEN 
GtSIREO BRANCH CUTS ANC OIRECTION ♦IMAGINARY CR “IMAGINARY 


FUNCTION S3RT1 (9ETA/Z) 

IMPLICIT OOU6LE PRECISION (A-H/O-r) 

CCMPUEX*U BeTA/AA/SORTI/CI 
COMMON /CONST/ CI/PI 
COMMON /CCNST1/SPEED/CMEGA 
COMMON /CCNST2 / AlFmA/CTCT 

AA * 1.-((CSPE£C/CMEGA)*8STA)«(CSPEE0/0FEGA)*8ETA)) 

o3 * 1. ♦ ALPHA*! ♦ OTOT 

S0RT1 »COSCRT (AA*2B-0TCT) 
ki TURN 

END ‘ 


* 
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n on fi o o n o c 


ORIGINAL PAGE 

c OF POOR QUALITY 

c 

c 

function sqrt2<beta) 

IMPLICIT DOUBLE PRECISION (A-H, 0 -Z> 

CCMPLE 1*16 9 ETA/AA,S 0 RT 2 ,CI 
COMMON /CONST/ CI*PI 

CCNNCN /C 0 NST 1 / SPEEC/OMEGA " 

COMMON /CON ST 2 / lt.PH»,CTCT 

AA * 1.~(((SPEcC/C ME GA)*9ETA)*(($ PEED/ OMEGA) •BETA) ) 
3 os t 2 * c:s;«t(aa) 
kc turn 

ENO 

c 

c 

c 

c 

SuSRCuTINi GALL(IsEcTAsG) 


C GALL c/AuLATcS T r- E j FUNCTION 


* 


IMPLICIT OCuELi PRECISION <A-M,0-Z> 

COMPLEX *16 3/ icTA/1322/32/CI/gT 
COMMON /CONST/ Cl, PI 
COMMON /CC.NST1/ < PE EC s OMEG A 
COMMON /C3NST2/ ALPHA, CTCT 
CAL. GIOALL C,?t TA/OI2) 

Pm; * 0ATAN2 CCIFAG(G32),CREAl<G32)) 
if (cM .;t. 0 . dC ) phi*pM-2*pl 

■» * < (cs»t*l;I2) )**(2.dO/3.dO>> * cd««p(<phi*2.dO/3.dC)»cl) 

Si * ; 

If (PHI .lE. -P 1/2. OC) THEN 

3*0* CC£*P(2.CC/3.00*PI*CI) 

ELSE 

0*5* C0EXPU.CC/3.00*PI*CI> 

ENC IF 
RETURN 

END " ~ 


SUBROUTINE GZALL1(I/6ETA,GZ,GZZ,EN) 

0ZALL1 CALCULATES ALL TmE PARTIAL DERIVATIVES • 

CF TME g FUNCTION. ~ * 


implicit coueue precisionca-h,o-z> 

COMPLEX *16 EcTA,GZ,C 2 Z,G/GB/GS 8 / SORTS', EN,C I sOUM 

COMMON /CONST/ CI*PI 

COMMON /C0NST1 /SPEEO^CMEGA 

COMMON /C0NST2/AL?MA,CT0T 

CALL GALL(Z, 8 ET 4 ,G) 

SI > 1.00 

A * 1.*ALPMA*ZmCT0T ' ~ ~ 

GZ » 2.00 • SI* ALPHA * SQRT1 (BETA,Z) / <3.*C0S0RT<C*A> ) 

C • 2.*ALFHA**3.0C*0T0T/(9*A**2.00) 

GZZ * C/<GZ*G>-.5*GZ**2.00/« 

T • (3. * OMEGA/ ( 2 . * ALPHA*SPEED) ) **(2.BC43.00> 
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n o r r» n r ' r. noon r. nr»n nnoo 


ORIGINAL PAGE IS 
OF POOR QUALITY 


tN * T*i 

RETURN 

END 


acsRDOTZNE hAII(Z/h2/*-21/H1,h11,IRE6) 


n-tU car a SU5P3UT IN* CS3AIR TC CAUCUUATE 1/3 ORDER 
► 4N*at. 'uNCTIONS FSCM AI5T FUNCTIONS. 


.mPlICIT COUsle PRECISION <A-«/0-Z) 

complex* is :,4:,6:/A:?/aiF,x,KS/Mi,H2,Nii#M2i 

uCM? l a X • 1 i HZ, ci 

uCMMuN /const/ CI/PI 

arc - ::m?lx<c.oO/-p:/ 6 .oc) 

.< * <ic. :•:)••(! .cc/6 .cc)«coexr<arc> 

* j * O.ONJGCA) 

C-UU CO: l-.zt-l, il, sI/A:P,5IP/!R?G> 

'i * *•< a;-c: *ai ) 

nC - *S # \AI*CI # :I) 

nil * -x«(a:»-c:*-::o> 

n21 a -<a«(AIP*CI»sIP) 

Pc Tj P *. 

END 


jUaSCUT INE C0eAIR<Z/AI/3I/AIP,SIP/IR«C> 


C-LuJUATE AIRY “UNCTIONS FCR CCMPUfX*1E ARGUMENT 

»:f . rANoecox op pat>ie«at:cal functions* asrafoniu ano stegun. 

ENTRY : 

CAUCucAT a IROjMaNTC) ANC AESCUUTE VAUUECZ) 

*P /!/ UT 6 

TncN USE cCS. 1C. 4. I THRU 1C. 4. 5 FOR A 1/ 3 1 , AIR, IIP 
C 1u uk Ir ASCII) UT PI/3 
C TnaN CAUCUUATE lETA(Z) 

C Uac e.S. 1C. 4.54/ 10.4.61/ 1C. 4. 63/ 10.4.66 FOR AI/SX/ AIP/BIF 

C cJ ciae CAUCUUATE ZeTA(-I) 

C USE aia. 10.4.60/ 1C. *..62/ 10.4.64/ 10.4.67 FOR AI/ 8 1/ AIF/BX P 

C cNOI r 

C cNDIP 

C cxIT 

U c N 0 


IMPUICIT 00UBU6 PRECISION (A-h/O-Z) 

COMMON /CONST/ CI/PI 

CCMPuEx»16 Z/AI/3I/AIP/EIP/ZETA/CZETA/Z14,SUN1/SUN2/SUN3*SUN4/ 
1 CETA?/ c ACT1/FACT2/SN/Ca/FTcRM/FPTERM/GTERF/GPTERM/F/FP/G/6P/Z3 
C0MFU;X»16 V 2 ET A / v Z ET AP,C I 
DIMENSION C C 2 1 )/0(21 ) 

CATA C1/C2/PIRT/P 14/. 35 50280539 00/. 25881 9403860/1. 77245385100/ 

♦ .735398163500/ 

PIRT*CSCRT(PI) 

DATA C/1.D0/.C69U444444444440C/ 

♦ .C 37 1 3343765 4 321 DC/ .C37993C591 278C000/ 

1 .C57 64 919C412669DC/.11 60 99 C 6 4C 25 51 00/ 
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ORIGINAL PM& IS 

♦ .291591 3992307400/. 877466969SC998C0/ OF POOR QUALITY 

2 3. 0794530 301731 00/12. 341 57 33 32 34SCO/ 

♦ 55. 62278536591400/278. 465C8C77759C0/ 

3 1533.16943 2012700/9207. 2 06 599725 8C0/ 

♦ 59892.51 356587500/419524*8751145300/ 

4 3U82S 7. 417866600/ 25198919. 871601C0/ 

♦ 214236036.9636600/1929375549.182300/ 

5 16335766937.86900/ 
j a r i c/1.00/ 

♦ -.09722222 2222221 CO/ -.C4 38 65 03C 864 19700/-.042462830789894C0/ 

1 -.0626621 6 3492031CC/-.1 241 C 58960 27270C/-.3082 53764901 0700/ 

2 -. 9204 79 r » 24 1 29 10C»- 3. 21 049358 4646 500/-1 2. 80729308073500/ 

3 -57. 5C32C3 51 39 110C /-28 7.023 2371092000/-1 576.357303337000/ 

4 -944s. 35 4e2 309 5 30C/ -all 3 5. 70666 384700/ -428952. 4004000400/ 

5 -321*536. 521 4CC60C/-256979C8. 363 9C 900/ -21 8293420.8321400/ 
e -1 9635 2 37 -i.99C90C/-1£ *43931088 .1C500/ 

4651 * 4*5(1) 

ZFUasZ.iC.O) GO TO 3 

;F<A3$to:“A&u)>.ie.i.c-i2.ANo.o#e*L<n.iT.o.co> co to 5 
4 S.l i CAT4N2<31*»G(I),0REAL(Z>> 

C TO A 

l ARGZ * c.:g 
«C TO * 

‘ ARGZ * D I 
I CONTINUE 

:?ca3$z.gt.6.co> cc tc io 

ir • ; = 1 

ASCENC1NG SERIES 
c*S • 1 Z • * • 2 0 1 C • A « 5 
CONTINUE 
Z3 * Z* *3 
FTcR* * 1.00 
FPTERP* Z*Z/2.00 
*T cRH * Z 

GFT ZZP* 1.00' ~ 

GCi* * 1.0-13*Aft3l 

F * FUR* 

F? * FPTc** ■*’ * 

G * GT ER* 

*P * GPTSR* 

KKKT * ICO "T "AC JUST KKKT TO INSURE CO*U*Gt*CI XP NECCSSARV 

CO 1 

13 « 3*1 

FTER* * FTcR«*Z3/C(I3-t.OC>«l3X 

FFTeRPsP?TERN*Z3/(I3*(l3^2.CO>) 
iTcSN * GTjRN/ZS/ClS/dS^I .00) 

GPTERP S GPTERH»Z3/{(I3”2. 00*13) " 

F * F+FTeRN 

FP » FP*FPT£RN 

C * /♦6TERM " " ’ - - 

iP * GP4CPTEBN 

ZF(C0A8S<CTERl»>.Lf .GLIN> GO TO 2 

1 CONTINUE * 

PRINT 6000/ 1 

C FORMAT { / * l«'2iU.5/* ERROR IN COB AIR/ NONCONVIRGENCE*) 

2 AI ‘ » C1*r-C2*G 

AZP • C1*8P-C2*GP 

81 • 1 .732C508C800* (Ct*MC2*G) 

8IP • 1.732C508C8D0*<CT*FFK2*CR3 

GO TO 9999 ’ . 
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ORIGINAL PA.GL is 

c asymptotic expansions PCS m large POOR QUALITY 

1C SIGN » 1 .30 * 

SUN1 • 0.00 
SUM 2 • 0.30 

3 U H 3 * C • 0 0 

SU H4 * C • 00 

If U3S(A*GZ) . GE • PI/3.00) GO TO 20 
C /i^C)/ Lc PI/3 

c eCS. 1C. 4.59,. 10. A. 61# 10.4.62/ 1C. 4. 66 

C : T A * CZETAUdSZ/ARGZ) 

-c ii : = i , i £ - 

k = 

C:TA ? - I £ T i • * « 

,w^1 * 5U'0*S:.iK*C(I)ZZSTAP 

■><.*2 * 3jM£/si/A«C(:)/2£T4P 

SUMS s 5wMi*C { I ) / ZETAP 
-o*w = 3„M4*i: ( ; ) / i£T a? 

11 L 3 n t ir u# 

;i, * a si -..is:: • ocmplxccosc6«ge/4.oO)/Sin<*rg:/4.co)) 

<-A»i 1 s .SCO* C3cXP(-Z£T»)/{»IRT*2l4) 
riCTi: * . 530* C3£xP(-ZETA W14/PIRT 

a: * *4CT1»SUN1 

-I* * -=ACT2*SLM2 

f*CTi * c:exp(:st*)/(p:rt*zu) 

r-CT2 s C3SXF(Z:T4)*Z14/PXRT 
si S ? AC T 1 • SUP I 

0 :? s P AC T2 » SUP4 

/ v T C 

/AnSCI)/ GT =>I /3 note CHANGE ABOVE 
£;s. 10.4. eC, 1C. 4.62# ic.4.64/ 1C. 4. 6? 


2C CONTINUE 
i r«g=3 

* oatan2(-3i*ao czj /-opeal cr) ) 

ZcTA * CZ=TA<AfiSZ/4RGZ) 

VZ=TA * 1.3C/ZcTA 
ILL * 10 
00 Z 1 :«1/LLL 

*2 * <x- 1>*2 

J * *2*1 “ “ 

VZcTAP * VIE TA**ic2 

SOMI * SUH1 ♦SIGN*C< J) *VZ5TAP 

SUM2 * SUH2*(SIGN*C ( j*i ) *vzetap*vzeta> 

SUH3 * 5UM3+($ZGN+C(J)*VZ£TAP) 

SUM4 * SUM4**CSIGN*C<JM)*VZETAF*VZETA> 

SIGN * -SIGN 
21 Continue 

m « A6SZ*«.Z500*OCMPLX(CCS<ARGZ/4.0C)/SIMARGZ/4.00)) 

FACT1 = 1 .30/ <PIRT*Z14) 

PACT2 * 214/PXRT 

SN » SINCZETAPPI/4.0C) 

CS * C0S(Z£TA*PI/4.0C) *' 

AI * PACT1»CSN*SU«1-CS*SCM2) 

AIP * -FACT2*(CS»SUM34SN*$UM4) 

az • facti*ccs*soni*sn*sum2) * • 

6IP * PACT2*(SN*SUM3-CS*SUM4) 
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°r*or*o oooooooo o n n n n 


c ORiG’MAL PASS IS, 

9999 continue ' OF POOR QUALITY 

RETURN 
CNO 


BEGIN OF FUNCTIONS USEC ABOVE 


function cieta ussz, *rgz> 

IMPLICIT DOUBLE PRECISION CA-N,0-Z> 

ccmplexob cieta 

ARC * AR S Z • 1 • SOC 

CIETA * (ABSZ**1. SuO^CC MPLR<CCS(ARG)/SIN(ARG) ) * .646666664 666 £ TOO 

ftiTuRN 
C N J 


cKC OF FUNCTIONS IS£0 ABOVE. 


iUBnOOTXNE *AIl 2< Z*H<, hi ,IREG> 


nAtLi uSeS SUfiaCuTISc CG 5 A I R TC CALCULATE 1/3 ORDER • 

h A NAc L s UNCT::nS PRCM AIRY FUNCTIONS. NCT THE * 

CcRIVmTIvES AS h d l L lCES. • 


IMPLICIT CQual: *<?ECISICN <A-H,C-Z> 

CChmCN /CONST/ CI/PI 

CC*=>l = x*1 6 :,a:*3I/AIF,3IP,X,XS*m1/M2*H11,H21 
CCMPLcX*1C ARG/CI 
CC*PLcX*l6 3 S T A 
AR G = DCMPLX(C.OC*-PI/6.CC) 

A * <12.CG)**<1 .00/4.00 *CCcXP(ARG) 

aS * OCCNJGU) 

CA w l CG5AlR(-:/AI/3I^AIP/eiP#IREG) 
n 1 * X*U:-CI*3I> 

«2 * KS*CAX+CZ*EX) 

RETURN 
SNu 

c 

c 

c 

C * - 

Subroutine mankel 

C 

C SUBROUTINE TO ACCURATELY CC THE HANK EL TRANSFORM CM THE SOUND 
C PRESSURE LEVEL. 

C F CNP) «GBAR<NP) NIST EE SAMPLEO AT NP POINTS WITH K-CN-1/ALP) 

C ALP REPRESENTS THE DISTANCE ABOVE THE REAL AXIS THE FUNCTION HILL 

C BE INTEGRATED. 

C NS IS A PARAMETER REPRESENTING AOOITION OP AN ANALYTICAL FUNCTION 

C TO FCNPJ " '• ’~ 

C M IS THE NUMBER CF TERNS USEO TO APPROXIMATE PCNP) TO INFINITY 

C ..........••**••••.*•*»«*•*••**•*•****••*•*•»•*•••**•* 

IMPLICIT DOUBLE PRECISION a-N,0-I> - 

COMMON /CONST/ CPPI/PI Vf 
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n o ft 


C 

C 


ORIGIN A L pfirte f 

0F p °°* qoalitv 


CCNMCN /."AIN/ p 

common /:ntsg/ns/me/*.i 

CONHCN/AFa2IN/ALF/Oel6eTA 
CGNPLEX.U ? (327*3) /CP/ CARG/SUN/PNP/CNPI/01 

At? * -ALP ! nacasiary to aatc* at tanborouflh'* dafioltion of alp 

N? x *1 
OeLK * DEL 3 ET A 

iteTSACT Tm£ ANALTTICAL FUNCTION IP NS > IfRO 
ACJuSTI.NG Tnj juETPACTIGN MULTIPLIER CP 
IF(N$.LE.J) GOTO 11 


C s » CCMPLX (C.CC/O.CC) 
IP (AtP.ti.C.C) 


cp « :?l:at(np)/:plcat(ns) 
CP x c = *m;) 


i m 3 IP 

IP UtP.NE.O.:) Tr;N 

CP * CMp;»CFtCAT ( NF ) • * ( 1 ) / (CPLCAT<NS)oAtP> 
e NO IP 


JtstS.CT T"; ANALYTICAL e LNCT ION IP NS>0 
CC K/Ixl/NP 

Cl * :CM?LX C e LCATC-1 ), (-ALP)) 

cars* cp t:AT<Ni>.(-:i)/;atCAT<NP) 

: < I ) * P CI)-C i »l 1 .CC-COEXP(CARG) ) 

1- CCNTlNUs 

11 C CNT INUE 

ip(Atp.*;. 3 .c) *<i ) s:cmplx<o.cO/C.oo) 

?N? s r(\ 0 ) 


. Z 12 /C* 2 /NP 

31 * 0 C-fL»< 0 PlCATC:- 1 ),<-AlP>> 

MU* = c)/(cos:*T(ci>> 

12 wwNTIiNw^ 

ipcalp.h*.:.c) *m*p n)/ccscRT((-CHPi)»(AtP)> 

c acc ii<«4 tc : x ; n : t t :* <*>o 

IP(*S.lT. 1 > i C Tj 20 
CC 1 5 / I si /NP 

31 * CC MB LAC 0 PtCAT<I- 1 )/(-AtP)) 

CP » 31 /C e tOAT (NP) 

CALL IcTA (NP, He/CP, SUM) 
r(I)s r ( 2 ) * F N P * S LM 

13 CCNTINUS 

23 CONTINUE 

C OC The rFT 

CAtt PCRKCNP/P/ 1 ) 

C ACC ALTERNATE TERMS TO GIVE NP /2 SAMPLES TRANSFORMED 

Cr x 35 LP*CPLCAT(NP)*(CCSC«T(-CMPI))/< 2 . 00 *PX> 
CO 25 / I * 2 / NP /2 

A s CeXP(CFL0AT(I-1)*(AtP)«2.C0*PI/0PL0AT(HP)) 
r ( I > x A«F(I)»ChPI»F(KP- 2 '* 2 )/A 
r(I)x f (I).CF/0*<.AT(CFLCAT(I-1)> 

23 CONTINUE 

RETURN 
ENO 


SUBROUTINE FCRK(lX,CX/SI 6 NI) 

C 

C A PAST PPT GIVEN it J.p. CtAERBOUT/ "FUNDAMENTALS OP GEOPHTSICAC • 
C LATA PROCESSING** • 
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noon o n r> n r- n o r or» 


ORIGINAL PAGE IS 
OF POOR QUALITY 

c*®*«» «...***•.****•• ••••«••••*•••••*••*•••**•*»«*••• 

IMPLICIT 00U8LE PRECISION <i-H,0-Z> 

COMMON /CONST/ CMPI,PI 

CC^LEX* U Cl CLX) »C ARG/CW/CTEMP/CI.OUMI /CMP I 
INTEGER SIGNI " 

J » 1 

.. * SC ■ OSQRTEI.OO/CFLOATiLX)) 

00 3C* I*T#LX 

IF(I.CT.J) GOTO 10 

CTcMP * CX(J)«SC 
CX ( JJ » CX(I)®SC 
CX<I) « CTcMP 
10 m*LX/2 

2u IP(J.Le.M) SC TO 3C 

J * J-M 

.1 * !*/ 2 

IFCM.Ge.l) SCTO 20 
iO J * J*M 

l * 1 

<.0 ISTrF*2»L 

CO 5C/M*1 /L 

CA«;=C«PI«®:*:6L;(S:ONI)«C3LE(<M-1))/OeLE(L) 

C-=CCEXPCCA?G) 

00 5CxI*M*LX, ISTe® 

CTiMPs Cm®CX(I«l> 

CXC*L)*CX(I)-CT:MP 
SO :XC)*CXC)*CTEMP 

L * I S T e P 

1 - (L .LT . LX) SCTO *C 
« C To i N 

•SC 

C 

c 


iJSSCLTINE I:TX(XP/M/*/SLM) 


S03SCLTINS TC »3C THE NECESSRRT TERMS TO 
jXPRESS INVER TIE LE FUNCTION TO INFINITY. 

■ Ill USE OOL o L i PRECISION. 

SUM® SUM OF 1 / (NP».5)*1/(<JPR}*.5) FOR J«1 TC 
INFINITY MINUS SOME CONSTANT MUCH IS 
INOcPENCeNT OF A. 


IMPLICIT OOUELE PRECISION CA-H,0-Z> / 

CCMFL2X®16 A/ SUN/C 2 

02 « OCMPLX (DFLC4TOO/C.CO) 

SUM ■ 1 .00/CMTA) 

SUM ■ 2.00*(OSORTCOFIOAT(M))-1.CO/COSQRT<SUM>I 
$ -C.5®CDSSRT<SUM)®(1.C*SLN*<1 .0/12.0-SUM*SUM/192.0J> 
CO 1C/J*1*M 

SUM « SUM+1 .OO/CCSQSTCJ+A) 

10 CONTINUE ' 

SUM ® SUM/OSQRTCCFLCATCNP)) 

RETURN 

EMO ‘ " 
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ORIGINAL F-3.E S3 
OF POC^. r^AS-'TY 


IjeiCLTINE 0RIHT4ll<CEU*Tl,N) Vi 

C 

C 

c • 

c » s*' 3 »ouTme to rriht the results or t>e 

c • r AHKcL TRANSFORM TO THE RILE R 0 RC 44 .CRT 

C • 


illicit cousie precision (a-h,o-z) 
c i H H c n / c c N j r / c I / s i 
cchhcs / , *a:n/ « 

CCHrLix-la e <2<7‘3),CI 
CO ICC/I*1/N/2 

'•»: * c.cc»f;« cpl:*tu-d / < crlo»tcm>*del»et* > 

Ip HAC .iT. C.O) This 

CiCIi.i * 2'.:c • :LC01C(CCA9S(R<I)>) 

^ £ = W> L V> L 1 » ( < i - ) 

= i*i,yz:) rac,rac 2 /:ecible 

i'*j l r 

1 iwNT »i* Jj 

r^C ? : * < - T ( : t , 5 G 1 5 . 7 ) 


sl5>:ct;nz p=g:ck. c :n: (anG/Rti ,rt2) 


c • T ms susrcutir# <J*ttrmin#$ there region changes 

c • t-e g2/2 function take place. The routine till 

c * :a : ? 1 1 »c by the earn pregrae anti eill return 

c • tne varieties » S 1 / R S 2 / RZ 1 * RZ 2 , RQ1 anti R02 ehlch 

- * 3r# th# values of beta ehare region changes occur 


illicit ccv?le precision ci-m,q-i) 

cc^mcn /const/ ci/Pi 

COMMON /CChSTI/ SP-EC/ChEGA 

CC*MC* /C0hSt2/ ALPnA,GTCT 

Cw^-hCN / C CN 3 T 3 / Z/IhcF^S 

CvMMCN /SEGICS/ CS1,R$2,SZ1,RZ2,R01/R02 

COMMCN /1 c 32IN/ hGHT, CELESTA 

CCMMCN / I h T r G / KS/hE/NDPTS 

CCNPtcX»l6 a, Cl 

Cl PENSION RSC3)/CZ(3)/R0(3) 

print • ✓ 'first rcot i$ # *RT1 

print */‘s»ccnd root # #RT2 

hLI.1 = * T 2 * 1.C100 

3LIH = RE4LUNT(.95C0*RT1/06L6ETA>) * 0EL6ETA 
5 N = H 3HT *DE L S ET A 

I X 0 ' ‘ ' * 

J = c 

K * C 

6 a OChPL X (BLIh/H) - 

1C CALL G32R6GICh(S/6/I/AMG#PS1) 
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ORIGINAL PAGE IS 

OF POOR QUALITY 


CALc G3t?cuICMZ/3/J«ANG*P!1) 

CALL G!2ft£GICN<0»B#R*ANG/P01) 

Aid) • REALd) 

HU) » REAL'S) 

ACU) « REAL (A) 

3 * 5 ♦ Zci.it! I 

If ((CREAtd) .Le. hl;<*) .ANC. (!♦.»♦« .LT. 9)> GOTO 10 

Ip C:*J*K .LT. iJPRIM • ,* REGION tOUNOART NOT OCTECTEOp* 

55 A $ 1 »« 5 ( 1 ) 

SS2««id) 

All *RI (1 ) 
ai2*RIU) 
aci*iC(i) 
iw2*dd) 

e nt * ;:.NT((N:N(ssi/a2<#aCt)-R0i)/06Le6TA) 

if irt . ; t . I ) t r «n 

p *InT *# *«su ■'*iri*ui n^foir cf point* in region 2 is**NT 
• no 1 1 

1 print *,'RS1,P$Z RSI, *32 
print •,'AZ1,*:2 # ,R!1, *Z2 

print */*3C1#5C2 # ,*C1,AC2 

pi c T u * N 
i NO 


Su3k;u t :n= G::RcG!0N<Z,B,L,ANG,P1) 


* SwC«:l:::ic TC CcTpRRlNE 332 FUNCTION ANO WFERE REGION CHANGES OCCUR * 
••••**•* 

implicit c:j = le ^rcisrcNca-^c-n 
ccmmcn /const/ c:,p: 

C 3 MPtcX* 1 o 2 , * 32 , Cl 
CALL 332ALLCZ,3,G32) 

?c * 0ATAN2<o:mag<G32),:R5AL(G32)) 

IF C?2 • GT • C.JJ) P2*P2-2*PI 
IF (L .EC. 0) THEN 

L * 1 

GOTO 1 
; NO IF 

IF C C r 1 .LT. ANG .ANC. P2 • SE • ANG) .OR. 

i CP 1 ,3c. AN 3 .ANC. P2 .LT. ANG)) L»L*1 

1 r 1 s ? 2 

RETURN 

£ NO ' 


SU3* CUT I Ne PE3I0N < BETA , IS EGICN) 


Thu suoroutin# ■ill dtttralnt «Mch region it • 

carrmt ly bting addrtsstd by tNt progrta * 


IMPLICIT OOUBLE FRECISI0N(A-H/0-Z> 

COMMON /C0NST3/ ZsZREF/S 

COMMON / REGION/ R $1 , R $ 2 , R 1 1 , RZ 2#R01 # *02 
COMMON /AFB2IN/ l-GHT, DELBET A 

complex^u bet* ’ ' 1 — : — ■ 

6 « 0tEAL<BETA)-0ELBITA/5.DO . , \ 


no 


ORIGINAL PAGE IS 
OF POOR QUALITY 


IF <2 .GT. S) THEN 


Iff EGION ■ 1 
(C .Iff. R02> ) Iff EGION * 6 
<e .LE * ff S 2 ) ) Iff EGION * 3 
(B .IE. ff 22 ) > Iff EG I ON » i 


IF (<d .GT. 
IF ((8 .GT. 
IF ((3 .GT. 
ELSE 

IF <(3 .GT. 
IF <<B .GT. 
IF UB .GT. 

SNC IF 
« cTuRN 
END 


ff 01 ) .AND • 
ffSD .ANC. 
ff 2 1 ) .AND. 


ROD .ANC. 
R21) .AND. 
RSI) .AND. 


(6 .16. ffC2>) 
(5 .LE. ff 2 2) ) 
(9 .LE. ff S 2 ) ) 


I REG ION » 2 
Iff EGION * 7 
I ff EGION > 8 
Iff EGION > 5 
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